TURB_LesMij.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 the mij stress appearing in the dynamic LES models.
!
! Description: We determine mij as appearing in the dynamic formulation:
!              M_{ij}=-test[bar(rho)](kappa \Delta)^2 |S(v)|S_{ij(v)
!                     +test[bar(rho)\Delta^{2}|S(u)|S_{ij}(u)]
!              in which v=test[bar(rho ui)]/test[bar(rho)] where test denotes
!              the test-filter with twice the filter-width. M_ij is computed
!              on the cell faces. S_ij(u) and Sij(test[u]) are already 
!              available, respectively, in mSij for mixture field and fSij for 
!              filtered field at cell faces. Moreover, we use the fact that we 
!              know 1/test[bar(rho)] (rhoFTb) and 1/[bar(rho)] (rhoF) at faces.
!
! Input: region = data of current region
!        ijk    = ijk-face is being treated
!
! Output: mij = resulting stress tensor in dynamic LES models
!
! Notes: none.
!
!******************************************************************************
!
! $Id: TURB_LesMij.F90,v 1.8 2008/12/06 08:44:41 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_lesmij( region,ijk )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
#ifdef RFLU
!  USE RFLU_ModInterpolation, ONLY: RFLU_InterpFaces2Faces
#endif
#ifdef RFLO
  USE modinterfaces, ONLY : rflo_getdimensdummynodes, rflo_getnodeoffset
  USE turb_modinterfaces, ONLY : turb_flolesunifiltff, turb_flolesgenfiltff

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

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

! ... loop variables
  INTEGER :: i, j, k, ijkn, ijkc0, ijkc1

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

  INTEGER      :: ibn,ien,idbeg,idend,tndel(diri:dirk)
  REAL(RFREAL) :: kappa2,rkappa2,two3rd
  REAL(RFREAL) :: delfac2,delfkap2,forterm
  REAL(RFREAL), POINTER :: fvol(:),fsij(:,:),fvar(:,:),ffvar(:,:)
  REAL(RFREAL), POINTER :: sij(:,:),mij(:,:)
  REAL(RFREAL), POINTER :: field(:,:),tfield(:,:)
  REAL(RFREAL), ALLOCATABLE :: modstrain(:)

#ifdef RFLO
  INTEGER      :: idnbeg,idnend,jdnbeg,jdnend,kdnbeg,kdnend
  INTEGER      :: ilev,inoff,ijnoff
#endif

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

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

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

! get indices and pointers --------------------------------------------------

#ifdef RFLO
  ilev = region%currLevel
  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
                                 jdnbeg,jdnend,kdnbeg,kdnend )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )

  fvar  => region%levels(ilev)%turb%fVar
  ffvar => region%levels(ilev)%turb%ffVar
  mij   => region%levels(ilev)%turb%mij
  ibn   =  indijk(idnbeg,jdnbeg,kdnbeg,inoff,ijnoff)
  ien   =  indijk(idnend,jdnend,kdnend,inoff,ijnoff)

  IF (ijk==diri) THEN
    fvol => region%levels(ilev)%turb%fvolI
    sij  => region%levels(ilev)%turb%mIsij
    fsij => region%levels(ilev)%turb%fISij
  ELSEIF (ijk==dirj) THEN
    fvol => region%levels(ilev)%turb%fvolJ
    sij  => region%levels(ilev)%turb%mJsij
    fsij => region%levels(ilev)%turb%fJSij
  ELSEIF (ijk==dirk) THEN
    fvol => region%levels(ilev)%turb%fvolK
    sij  => region%levels(ilev)%turb%mKsij
    fsij => region%levels(ilev)%turb%fKSij
  ENDIF
#endif
#ifdef RFLU
  fvar  => region%turb%fVar
  ffvar => region%turb%ffVar
  mij   => region%turb%mij
  fvol  => region%turb%fvolI
  sij   => region%turb%mIsij
  fsij  => region%turb%fISij
  ibn   =  1
  ien   =  region%grid%nFaces
#endif

! allocate work arrays used within this scope

  ALLOCATE( modstrain(ibn:ien) )
  ALLOCATE( field(e11:e33,ibn:ien),tfield(e11:e33,ibn:ien) )

! test filter width is twice bar filter width

  tndel(diri) = 2*region%turbInput%filterWidth(diri)
#ifdef RFLO
  tndel(dirj) = 2*region%turbInput%filterWidth(dirj)
  tndel(dirk) = 2*region%turbInput%filterWidth(dirk)
#endif

! filterwidth on the fg-level (test-bar) is sqrt(5) times larger than 
! on the f-level (bar), i.e. kappa=sqrt(5)

  kappa2 = 5._rfreal
  rkappa2= 1._rfreal/kappa2
  two3rd = 2._rfreal/3._rfreal

! determine the bar filter width delta which appears in the dynamic model; 
! this factor in square is defined as (kappa*delfac)^2 and stored in delFKap2 

  delfac2  = region%turbInput%delFac2
  delfkap2 = delfac2*kappa2

! first calculate test[bar(rho)](\kappa \Delta)^2|S(v)|Sij(v) at cell faces; 
! we obtain Sij(v) from array fSij; 

  DO ijkn=ibn,ien

! - generate the 6-th component of fSij using the tracelessness of S_{ij} 

    fsij(e33,ijkn)=-(fsij(e11,ijkn)+fsij(e22,ijkn))

! - get modulus of Sij(v): |S(v)|
    modstrain(ijkn)=sqrt(fsij(e11,ijkn)*fsij(e11,ijkn) &
                        +fsij(e22,ijkn)*fsij(e22,ijkn) &
                        +fsij(e11,ijkn)*fsij(e22,ijkn) &
                        +fsij(e12,ijkn)*fsij(e12,ijkn) &
                        +fsij(e13,ijkn)*fsij(e13,ijkn) &
                        +fsij(e23,ijkn)*fsij(e23,ijkn))

! - in turn, (kappa*delta)^2 in fg-level is stored in tField

    tfield(e11,ijkn) = delfkap2*fvol(ijkn)**two3rd

! - fg-level contribution to M_{ij} can now be computed since
!   we already have test(bar[rho]) available at cell faces,
!   in fact we have 1/test(bar[rho]), so we devide instead of multiply;
!   install the first term of mij; the space of fSij in used in the process

    forterm       = -tfield(e11,ijkn)/ffvar(cv_turb_dens,ijkn)* &
                     modstrain(ijkn)

    mij(e11,ijkn) = forterm*fsij(e11,ijkn)
    mij(e12,ijkn) = forterm*fsij(e12,ijkn)
    mij(e13,ijkn) = forterm*fsij(e13,ijkn)
    mij(e22,ijkn) = forterm*fsij(e22,ijkn)
    mij(e23,ijkn) = forterm*fsij(e23,ijkn)
    mij(e33,ijkn) = forterm*fsij(e33,ijkn)

! - as final step we determine contribution of test-filtered f-level term;
!   note we have bar[rho] available at cell faces bRho stored in fVar; 
!   in fact we have 1/bRho in fVar(CV_TURB_DENS,:), so we devide by bRho
!   instead of multiply and compute modulus of Sij: |S(u)|

    modstrain(ijkn) = sqrt(sij(e11,ijkn)*sij(e11,ijkn) &
                          +sij(e22,ijkn)*sij(e22,ijkn) &
                          +sij(e11,ijkn)*sij(e22,ijkn) &
                          +sij(e12,ijkn)*sij(e12,ijkn) &
                          +sij(e13,ijkn)*sij(e13,ijkn) &
                          +sij(e23,ijkn)*sij(e23,ijkn))

! - finish computation of mij; filterwidth^2 at f-level is tField/kappa2;
!   build terms need to be filtered

    forterm = tfield(e11,ijkn)/fvar(cv_turb_dens,ijkn)*rkappa2* &
              modstrain(ijkn)

    field(e11,ijkn) =  forterm*sij(e11,ijkn)
    field(e12,ijkn) =  forterm*sij(e12,ijkn)
    field(e13,ijkn) =  forterm*sij(e13,ijkn)
    field(e22,ijkn) =  forterm*sij(e22,ijkn)
    field(e23,ijkn) =  forterm*sij(e23,ijkn)
    field(e33,ijkn) = -forterm*(sij(e11,ijkn)+sij(e22,ijkn))
  ENDDO

! perform test filtering:

  idbeg = e11
  idend = e33
#ifdef RFLO
  IF (region%turbInput%filterType == filtype_uniform) then
    CALL turb_flolesunifiltff( region,ijk,tndel,idbeg,idend,field,tfield )
  ELSE
    CALL turb_flolesgenfiltff( region,ijk,tndel,idbeg,idend,field,tfield )
  ENDIF
#endif
#ifdef RFLU
!  CALL RFLU_InterpFaces2Faces( region,field,tField )
  tfield(idbeg:idend,:) = field(idbeg:idend,:)
#endif

! and complete components of mij:

  DO ijkn=ibn,ien
    mij(e11,ijkn) = mij(e11,ijkn) + tfield(e11,ijkn)
    mij(e12,ijkn) = mij(e12,ijkn) + tfield(e12,ijkn)
    mij(e13,ijkn) = mij(e13,ijkn) + tfield(e13,ijkn)
    mij(e22,ijkn) = mij(e22,ijkn) + tfield(e22,ijkn)
    mij(e23,ijkn) = mij(e23,ijkn) + tfield(e23,ijkn)
    mij(e33,ijkn) = mij(e33,ijkn) + tfield(e33,ijkn)
  ENDDO

! deallocate temporary arrays

  DEALLOCATE( modstrain,field,tfield )

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

  CALL deregisterfunction( global )

END SUBROUTINE turb_lesmij

!******************************************************************************
!
! RCS Revision history:
!
! $Log: TURB_LesMij.F90,v $
! Revision 1.8  2008/12/06 08:44:41  mtcampbe
! Updated license.
!
! Revision 1.7  2008/11/19 22:17:54  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.6  2005/12/31 23:29:02  wasistho
! copy in lack of filter in rocflu temporarily
!
! Revision 1.5  2005/12/28 03:24:39  wasistho
! simply copy variables in waiting for Rocflu's FF-interpolation
!
! Revision 1.4  2004/05/28 02:01:12  wasistho
! update unstructured grid LES
!
! Revision 1.3  2004/03/19 02:50:37  wasistho
! prepared for RFLU
!
! Revision 1.2  2004/03/12 02:55:35  wasistho
! changed rocturb routine names
!
! Revision 1.1  2004/03/05 04:37:00  wasistho
! changed nomenclature
!
! 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
!
!
!******************************************************************************