RFLU_ComputeEnerDissOLES.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 kinetic energy and dissipation for use in optimal LES 
!   weight computation.
!
! Description: None.
!
! Input:
!   region      Region for which energy and dissipation is to be computed
!
! Output: None.
!
! Notes: 
!   1. At present, works only for single region
!
!******************************************************************************
!
! $Id: RFLU_ComputeEnerDissOLES.F90,v 1.5 2008/12/06 08:44:29 mtcampbe Exp $
!
! Copyright: (c) 2002 by the University of Illinois
!
!******************************************************************************

SUBROUTINE rflu_computeenerdissoles(region)

  USE moddatatypes
  USE moderror
  USE modglobal, ONLY: t_global
  USE moddatastruct, ONLY: t_region
  USE modparameters
  USE modtools, ONLY: makenonzero

  USE rflu_modoles
   
  IMPLICIT NONE

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

! ... loop variables
  INTEGER :: ic,ifc,iv1,iv2

! ... local variables
  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER :: a,c1g,c2g,d,errorflag,hloc1,ic1l,ic2l,ifcp,j,l,ncells,nfaces,vloc1
  INTEGER, DIMENSION(:), POINTER :: f2fpoles
  INTEGER, DIMENSION(:,:), POINTER :: fsoles
  REAL(RFREAL) :: avgfac,corr,eneroles,numer,sum1,sum2,sum3,term,uavg,var, & 
                  vavg,wavg
  REAL(RFREAL), DIMENSION(:), POINTER :: vol  
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv
  REAL(RFREAL), DIMENSION(:,:,:), ALLOCATABLE :: int2oles 
  REAL(RFREAL), DIMENSION(:,:,:), POINTER :: int20oles,int21oles   
  TYPE(t_global), POINTER :: global

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

  rcsidentstring = '$RCSfile: RFLU_ComputeEnerDissOLES.F90,v $ $Revision: 1.5 $'

  global => region%global

  CALL registerfunction(global,'RFLU_ComputeEnerDissOLES',&
  'RFLU_ComputeEnerDissOLES.F90')

! ******************************************************************************
! Start, check that have primitive variables
! ******************************************************************************

  IF ( region%mixt%cvState == cv_mixt_state_cons ) THEN 
    CALL errorstop(global,err_cv_state_invalid,__line__)
  END IF ! region

! ******************************************************************************
! Set pointers and variables
! ******************************************************************************

  vol => region%grid%vol
  cv  => region%mixt%cv

  eneroles = global%enerOLES

  fsoles   => region%grid%fsOLES
  f2fpoles => region%grid%f2fpOLES

  int20oles => region%grid%int20OLES
  int21oles => region%grid%int21OLES

  ncells = region%grid%nCells
  nfaces = region%grid%nFaces  
  avgfac = 1.0_rfreal/REAL(ncells,kind=rfreal)

! ******************************************************************************
! Compute kinetic energy (average over cells)
! ******************************************************************************

  eneroles = 0.0_rfreal

  DO ic = 1,ncells
    eneroles = eneroles + cv(cv_mixt_xvel,ic)*cv(cv_mixt_xvel,ic) & 
                        + cv(cv_mixt_yvel,ic)*cv(cv_mixt_yvel,ic) &
                        + cv(cv_mixt_zvel,ic)*cv(cv_mixt_zvel,ic) 
  END DO ! ic

  global%enerOLES = avgfac*0.5_rfreal*eneroles  

! DEBUG
  uavg = 0.0_rfreal
  vavg = 0.0_rfreal
  wavg = 0.0_rfreal

  DO ic = 1,ncells
    uavg = uavg + cv(cv_mixt_xvel,ic)*cv(cv_mixt_xvel,ic)
    vavg = vavg + cv(cv_mixt_yvel,ic)*cv(cv_mixt_yvel,ic)
    wavg = wavg + cv(cv_mixt_zvel,ic)*cv(cv_mixt_zvel,ic)        
  END DO ! ic 
  
  uavg = avgfac*uavg
  vavg = avgfac*vavg
  wavg = avgfac*wavg
  
  global%uVarOLES = uavg
  global%vVarOLES = vavg
  global%wVarOLES = wavg       
! END DEBUG

! ******************************************************************************
! Compute dissipation rate
! ******************************************************************************

  ncells = SIZE(region%grid%fsOLES,1)
  avgfac = 3.0_rfreal/REAL(nfaces,kind=rfreal)

  ALLOCATE(int2oles(3,3*ncells,3*ncells),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'int2OLES')
  END IF ! global%error

  int2oles(:,:,:) = 0.0_rfreal

! ==============================================================================
! Compute integral through loop over faces
! ==============================================================================

  DO ifc = 1,nfaces
    ifcp = f2fpoles(ifc) ! get prototype face

! ------------------------------------------------------------------------------
!   Loop over cells
! ------------------------------------------------------------------------------
 
    DO ic1l = 1,ncells
      c1g = fsoles(ic1l,ifc)
      d = ic1l
    
      DO ic2l = 1,ncells
        c2g = fsoles(ic2l,ifc)      
        a = ic2l
                         
        term = vol(c1g)*vol(c2g)

! ----- Loop over velocity components ------------------------------------------

        DO iv1 = 1,3
          l = iv1

          DO iv2 = 1,3
            j = iv2              

! --------- Compute correlation              

            corr = term*cv(iv1+1,c1g)*cv(iv2+1,c2g)

! --------- Determine storage location for correlation              

            vloc1 = rflu_geti1posoles(l,d)              
            hloc1 = rflu_getlposoles(j,a)

! --------- Store correlation              

            int2oles(ifcp,vloc1,hloc1) = int2oles(ifcp,vloc1,hloc1) + corr                                            

          END DO ! iv2
        END DO ! iv1
        
      END DO ! icl2
    END DO ! icl1

  END DO ! ifcp  

! ==============================================================================
! Normalize and average integrals
! ==============================================================================
                                      
  int2oles(:,:,:) = avgfac*int2oles(:,:,:)
   
! ==============================================================================
! Attempt to extract dissipation from I20 and I21 
! ==============================================================================   
   
  sum1 = 0.0_rfreal
  sum2 = 0.0_rfreal
  sum3 = 0.0_rfreal 
   
  DO ifcp = 1,3
 
    DO ic1l = 1,ncells
      d = ic1l
    
      DO ic2l = 1,ncells
        a = ic2l
      
! 1x1x2 stencil      
        IF ( ic1l /= ic2l ) THEN
! 1x1x4 stencil, innermost cells - should give the same as 1x1x2 
!        IF ( (ic1l == 2 .OR. ic1l == 3) .AND. (ic2l == 3 .OR. ic2l == 2) .AND. & 
!             (ic1l /= ic2l) ) THEN
! 1x1x4 stencil, outermost cells
!        IF ( (ic1l == 1 .OR. ic1l == 4) .AND. (ic2l == 4 .OR. ic2l == 1) .AND. & 
!             (ic1l /= ic2l) ) THEN
! 1x1x6 stencil, innermost cells  
!        IF ( (ic1l == 3 .OR. ic1l == 4) .AND. (ic2l == 4 .OR. ic2l == 3) .AND. & 
!             (ic1l /= ic2l) ) THEN 
! 1x1x6 stencil, 2nd-innermost cells  
!        IF ( (ic1l == 2 .OR. ic1l == 5) .AND. (ic2l == 5 .OR. ic2l == 2) .AND. & 
!             (ic1l /= ic2l) ) THEN                                     
! 1x1x6 stencil, outermost cells  
!        IF ( (ic1l == 1 .OR. ic1l == 6) .AND. (ic2l == 6 .OR. ic2l == 1) .AND. & 
!             (ic1l /= ic2l) ) THEN
! 1x1x8 stencil, innermost cells  
!        IF ( (ic1l == 4 .OR. ic1l == 5) .AND. (ic2l == 5 .OR. ic2l == 4) .AND. & 
!             (ic1l /= ic2l) ) THEN  
! 1x1x8 stencil, 2nd-innermost cells
!        IF ( (ic1l == 3 .OR. ic1l == 6) .AND. (ic2l == 6 .OR. ic2l == 3) .AND. & 
!             (ic1l /= ic2l) ) THEN
! 1x1x8 stencil, 3rd-innermost cells             
!        IF ( (ic1l == 2 .OR. ic1l == 7) .AND. (ic2l == 7 .OR. ic2l == 2) .AND. & 
!             (ic1l /= ic2l) ) THEN                                          
! 1x1x8 stencil, outermost cells             
!        IF ( (ic1l == 1 .OR. ic1l == 8) .AND. (ic2l == 8 .OR. ic2l == 1) .AND. & 
!             (ic1l /= ic2l) ) THEN     
          DO iv1 = 1,3
            l = iv1

            iv2 = iv1        
            j   = iv2

            vloc1 = rflu_geti1posoles(l,d)              
            hloc1 = rflu_getlposoles(j,a)

            sum1 = sum1 + int2oles(ifcp,vloc1,hloc1)
            sum2 = sum2 + int20oles(ifcp,vloc1,hloc1)
            sum3 = sum3 + int21oles(ifcp,vloc1,hloc1)            
          END DO ! iv1          
        END IF ! ic1l
        
      END DO ! ic2l  
    END DO ! ic1l 

  END DO ! ifcp
   
! NOTE: Assume that rhoOLES has a constant value throughout solution region, so
! use rhoOLES(1). This will need to be modified if the grids become non-uniform    

  var   = 2.0_rfreal*global%enerOLES/3.0_rfreal
  numer = sum1/region%grid%deltaOLES**6 - var*sum2
  global%dissOLES = (max(numer/sum3,tiny(1.0_rfreal)))**(3.0_rfreal/2.0_rfreal)/ & 
                    region%grid%rhoOLES(1)

  WRITE(*,'(3(2X,E15.8))') sum1/region%grid%deltaOLES**6,var*sum2,sum3
  WRITE(71,'(4(2X,E15.8))') global%currentTime,sum1/region%grid%deltaOLES**6, & 
                            var*sum2,sum3

! DEBUG
  WRITE(*,'(5(2X,E15.8))') global%enerOLES,global%dissOLES, & 
                           uavg/(makenonzero(3.0_rfreal*var)), & 
                           vavg/(makenonzero(3.0_rfreal*var)), & 
                           wavg/(makenonzero(3.0_rfreal*var))                                       
! END DEBUG


  DEALLOCATE(int2oles,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'int2OLES')
  END IF ! global%error

! ******************************************************************************
! End
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE rflu_computeenerdissoles

!******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ComputeEnerDissOLES.F90,v $
! Revision 1.5  2008/12/06 08:44:29  mtcampbe
! Updated license.
!
! Revision 1.4  2008/11/19 22:17:42  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.3  2004/01/22 16:04:33  haselbac
! Changed declaration to eliminate warning on ALC
!
! Revision 1.2  2002/10/08 15:49:29  haselbac
! {IO}STAT=global%error replaced by {IO}STAT=errorFlag - SGI problem
!
! Revision 1.1  2002/09/09 15:39:57  haselbac
! Initial revision
!
!******************************************************************************