RFLO_ModGridMetrics.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: Suite for grid metrics related routines.
!
! Description: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: RFLO_ModGridMetrics.F90,v 1.7 2008/12/06 08:44:16 mtcampbe Exp $
!
! Copyright: (c) 2004 by the University of Illinois
!
! ******************************************************************************

MODULE rflo_modgridmetrics

  USE modglobal, ONLY    : t_global 
  USE moddatastruct, ONLY: t_region  
  USE modgrid, ONLY      : t_grid
  USE modbndpatch, ONLY  : t_patch
  USE modparameters
  USE moddatatypes
  USE moderror
  USE modmpi
  
  IMPLICIT NONE

  PRIVATE
  PUBLIC :: rflo_arclengthpatch, &
            rflo_calcgridmetrics, &
            rflo_gridqualityglobal

! private :
    
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
   
  CHARACTER(CHRLEN) :: RCSIdentString = & 
    '$RCSfile: RFLO_ModGridMetrics.F90,v $ $Revision: 1.7 $'        
             
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS


!******************************************************************************
!
! Purpose: calculate approximate arclengths for every grid line
!          between two oposite patch boundaries (on the finest
!          grid only).
!
! Description: none.
!
! Input: region = grid dimensions
!        patch  = current patch
!        xyzRef = reference coordinates.
!
! Output: arcLen1 = arclength in first  coordinate direction of patch
!         arcLen2 = arclength in second coordinate direction of patch
!
! Notes: none.
!
!******************************************************************************

SUBROUTINE rflo_arclengthpatch( region,patch,xyzRef )

  USE modinterfaces, ONLY : rflo_getnodeoffset, rflo_getpatchindicesnodes

  IMPLICIT NONE

#include "Indexing.h"

! ... parameters

  TYPE(t_region)         :: region
  TYPE(t_patch), POINTER :: patch
  REAL(RFREAL), POINTER  :: xyzref(:,:)

! ... loop variables
  INTEGER :: l1, l2

! ... local variables
  INTEGER :: ilev, lbound, ibeg, iend, jbeg, jend, kbeg, kend
  INTEGER :: l1b, l1e, l2b, l2e, lc, ijkn, ijkn1, ijkn2, inoff, ijnoff
  INTEGER :: k1, k2, switch(6,5)

  REAL(RFREAL), POINTER :: arclen1(:), arclen2(:)

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

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

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

  lbound = patch%lbound
  ilev   = 1

  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
  CALL rflo_getpatchindicesnodes( region,patch,ilev, &
                                  ibeg,iend,jbeg,jend,kbeg,kend )

  arclen1 => patch%arclen1
  arclen2 => patch%arclen2

! set boundary switch ---------------------------------------------------------
! switch(:,1-2) = first/last index in l1-direction
! switch(:,3-4) = first/last index in l2-direction
! switch(:,  5) = constant index

  switch(1,:) = (/jbeg, jend, kbeg, kend, ibeg/)
  switch(2,:) = (/jbeg, jend, kbeg, kend, iend/)
  switch(3,:) = (/kbeg, kend, ibeg, iend, jbeg/)
  switch(4,:) = (/kbeg, kend, ibeg, iend, jend/)
  switch(5,:) = (/ibeg, iend, jbeg, jend, kbeg/)
  switch(6,:) = (/ibeg, iend, jbeg, jend, kend/)

  l1b = switch(lbound,1)
  l1e = switch(lbound,2)
  l2b = switch(lbound,3)
  l2e = switch(lbound,4)
  lc  = switch(lbound,5)

! compute arclengths of current patch ------------------------------------------

  arclen1(:) = 0._rfreal
  arclen2(:) = 0._rfreal

  DO l2=l2b+1,l2e
    k2 = l2-l2b+1
    DO l1=l1b+1,l1e
      k1 = l1-l1b+1
      IF (lbound==1 .OR. lbound==2) THEN
        ijkn      = indijk(lc,l1    ,l2    ,inoff,ijnoff)
        ijkn1     = indijk(lc,l1-1  ,l2    ,inoff,ijnoff)
        ijkn2     = indijk(lc,l1    ,l2-1  ,inoff,ijnoff)
      ELSEIF (lbound==3 .OR. lbound==4) THEN
        ijkn      = indijk(l2    ,lc,l1    ,inoff,ijnoff)
        ijkn1     = indijk(l2    ,lc,l1-1  ,inoff,ijnoff)
        ijkn2     = indijk(l2-1  ,lc,l1    ,inoff,ijnoff)
      ELSEIF (lbound==5 .OR. lbound==6) THEN
        ijkn      = indijk(l1    ,l2    ,lc,inoff,ijnoff)
        ijkn1     = indijk(l1-1  ,l2    ,lc,inoff,ijnoff)
        ijkn2     = indijk(l1    ,l2-1  ,lc,inoff,ijnoff)
      ENDIF
      arclen1(k2) = arclen1(k2) + &
            sqrt((xyzref(xcoord,ijkn)-xyzref(xcoord,ijkn1))**2 + &
                 (xyzref(ycoord,ijkn)-xyzref(ycoord,ijkn1))**2 + &
                 (xyzref(zcoord,ijkn)-xyzref(zcoord,ijkn1))**2)
      arclen2(k1) = arclen2(k1) + &
            sqrt((xyzref(xcoord,ijkn)-xyzref(xcoord,ijkn2))**2 + &
                 (xyzref(ycoord,ijkn)-xyzref(ycoord,ijkn2))**2 + &
                 (xyzref(zcoord,ijkn)-xyzref(zcoord,ijkn2))**2)
    ENDDO ! l1
  ENDDO   ! l2

  DO l2=l2b,l2b
    DO l1=l1b+1,l1e
      IF (lbound==1 .OR. lbound==2) THEN
        ijkn      = indijk(lc,l1    ,l2    ,inoff,ijnoff)
        ijkn1     = indijk(lc,l1-1  ,l2    ,inoff,ijnoff)
      ELSEIF (lbound==3 .OR. lbound==4) THEN
        ijkn      = indijk(l2    ,lc,l1    ,inoff,ijnoff)
        ijkn1     = indijk(l2    ,lc,l1-1  ,inoff,ijnoff)
      ELSEIF (lbound==5 .OR. lbound==6) THEN
        ijkn      = indijk(l1    ,l2    ,lc,inoff,ijnoff)
        ijkn1     = indijk(l1-1  ,l2    ,lc,inoff,ijnoff)
      ENDIF
      arclen1(1) = arclen1(1) + &
            sqrt((xyzref(xcoord,ijkn)-xyzref(xcoord,ijkn1))**2 + &
                 (xyzref(ycoord,ijkn)-xyzref(ycoord,ijkn1))**2 + &
                 (xyzref(zcoord,ijkn)-xyzref(zcoord,ijkn1))**2)
    ENDDO ! l1
  ENDDO   ! l2

  DO l1=l1b,l1b
    DO l2=l2b+1,l2e
      IF (lbound==1 .OR. lbound==2) THEN
        ijkn      = indijk(lc,l1    ,l2    ,inoff,ijnoff)
        ijkn2     = indijk(lc,l1    ,l2-1  ,inoff,ijnoff)
      ELSEIF (lbound==3 .OR. lbound==4) THEN
        ijkn      = indijk(l2    ,lc,l1    ,inoff,ijnoff)
        ijkn2     = indijk(l2-1  ,lc,l1    ,inoff,ijnoff)
      ELSEIF (lbound==5 .OR. lbound==6) THEN
        ijkn      = indijk(l1    ,l2    ,lc,inoff,ijnoff)
        ijkn2     = indijk(l1    ,l2-1  ,lc,inoff,ijnoff)
      ENDIF
      arclen2(1) = arclen2(1) + &
            sqrt((xyzref(xcoord,ijkn)-xyzref(xcoord,ijkn2))**2 + &
                 (xyzref(ycoord,ijkn)-xyzref(ycoord,ijkn2))**2 + &
                 (xyzref(zcoord,ijkn)-xyzref(zcoord,ijkn2))**2)
    ENDDO ! l2
  ENDDO   ! l1

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

  CALL deregisterfunction( region%global )

END SUBROUTINE rflo_arclengthpatch


!******************************************************************************
!
! Purpose: calculate face vectors, volumes, cell centroids (optionally),
!          and cell-to-face averaging coefficients
!
! Description: none.
!
! Input: regions%grid = dimensions, grid coordinates.
!
! Output: regions%grid = face vectors (si,sj,sk), volumes (vol),
!                        cell centroids (cofg), avg coeffs. (avgCoI,J,K).
!
! Notes: none.
!
!******************************************************************************

SUBROUTINE rflo_calcgridmetrics( regions )

  USE modinterfaces, ONLY : rflo_calcfacevectors, rflo_calccontrolvolumes, &
             rflo_calccellcentroids, rflo_calcfacecentroids, &
             rflo_initavgcoeffs, rflo_c2favgcoeffs, rflo_c2eavgcoeffs, &
             rflo_checkmetrics
  IMPLICIT NONE

! ... parameters
  TYPE (t_region), POINTER :: regions(:)

! ... loop variables
  INTEGER :: ireg

! ... local variables
  REAL(RFREAL) :: skewmin
  TYPE (t_global), POINTER :: global

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

  global => regions(1)%global
  CALL registerfunction( global,'RFLO_CalcGridMetrics',&
       'RFLO_ModGridMetrics.F90' )

! initialization some parameters ----------------------------------------------

  global%skewness = 1._rfreal

! loop over all regions -------------------------------------------------------

  DO ireg=1,regions(1)%global%nRegions
    IF (regions(ireg)%procid==regions(ireg)%global%myProcid & ! reg active and
        .AND. regions(ireg)%active==active) THEN              ! on my processor

      CALL rflo_calcfacevectors( regions(ireg) )

      CALL rflo_calccontrolvolumes( regions(ireg) )

      CALL rflo_calccellcentroids( regions(ireg) )

      CALL rflo_calcfacecentroids( regions(ireg) )

      CALL rflo_initavgcoeffs( regions(ireg) )

      IF (regions(ireg)%mixtInput%faceEdgeAvg==fe_avg_linear) &
        CALL rflo_c2favgcoeffs( regions(ireg) )

      CALL rflo_c2eavgcoeffs( regions(ireg) )

! --- check metrics

      CALL rflo_checkmetrics( ireg,regions(ireg) )

    ENDIF     ! region on this processor and active
  ENDDO       ! iReg

! global grid quality measure -------------------------------------------------

  CALL rflo_gridqualityglobal( regions )

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

  CALL deregisterfunction( global )

END SUBROUTINE rflo_calcgridmetrics


!******************************************************************************
!
! Purpose: global reduction of grid quality measures
!
! Description: none.
!
! Input: regions%grid, regions%global = grid quality data
!
! Output: global%skewness, etc = global skewness, etc
!
! Notes: none.
!
!******************************************************************************

SUBROUTINE rflo_gridqualityglobal( regions )

  IMPLICIT NONE

! ... parameters
  TYPE (t_region), POINTER :: regions(:)

! ... loop variables
  INTEGER :: ireg

! ... local variables
  REAL(RFREAL) :: skewmin, volmin
  TYPE (t_global), POINTER :: global

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

  global => regions(1)%global
  CALL registerfunction( global,'RFLO_GridQualityGlobal',&
       'RFLO_ModGridMetrics.F90' )

! global skewness and minVol --------------------------------------------------

  skewmin = global%skewness
  volmin  = global%minVol

  DO ireg=1,global%nRegions
    IF (regions(ireg)%procid==global%myProcid .AND. &     ! region active and
        regions(ireg)%active==active) THEN                ! on my processor
      skewmin = min( skewmin,regions(ireg)%levels(1)%grid%skewness )
      volmin  = min( volmin, regions(ireg)%levels(1)%grid%minVol   )
    ENDIF   ! active
  ENDDO     ! iReg

  global%skewness = skewmin
  global%minVol   = volmin

#ifdef MPI
  CALL mpi_allreduce( global%skewness,skewmin,1,mpi_rfreal,mpi_min, &
                      global%mpiComm,global%mpierr )
  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,&
       __line__ )
  global%skewness = skewmin

  CALL mpi_allreduce( global%minVol,volmin,1,mpi_rfreal,mpi_min, &
                      global%mpiComm,global%mpierr )
  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,&
       __line__ )
  global%minVol = volmin
#endif

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

  CALL deregisterfunction( global )

END SUBROUTINE rflo_gridqualityglobal


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

END MODULE rflo_modgridmetrics


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLO_ModGridMetrics.F90,v $
! Revision 1.7  2008/12/06 08:44:16  mtcampbe
! Updated license.
!
! Revision 1.6  2008/11/19 22:17:27  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.5  2006/03/24 00:53:51  wasistho
! fixed loop indexing in patch%arclen1,2
!
! Revision 1.4  2006/03/18 11:01:23  wasistho
! moved some routines to ModGridRegionShape
!
! Revision 1.3  2006/03/18 08:17:01  wasistho
! added arcLengthPatch
!
! Revision 1.2  2006/03/15 06:37:46  wasistho
! added region and global skewness
!
! Revision 1.1  2006/03/04 04:36:41  wasistho
! initial import RFLO_ModGridMetrics
!
!
!
! ******************************************************************************