RFLU_ModDimensionality.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 of routines related to two dimensional computations.
!
! Description: None.
!
! Notes: None. 
!
! ******************************************************************************
!
! $Id: RFLU_ModDimensionality.F90,v 1.9 2008/12/06 08:44:21 mtcampbe Exp $
!
! Copyright: (c) 2005 by the University of Illinois
!
! ******************************************************************************

MODULE rflu_moddimensionality

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

  IMPLICIT NONE
   
  PRIVATE
  PUBLIC :: rflu_123d_checkgeometrywrapper, & 
            rflu_123d_checktopology
      
! ******************************************************************************
! Declarations and definitions
! ******************************************************************************  
     
  CHARACTER(CHRLEN), PRIVATE :: & 
    RCSIdentString = '$RCSfile: RFLU_ModDimensionality.F90,v $ $Revision: 1.9 $'        
       
! ******************************************************************************
! Routines
! ******************************************************************************

  CONTAINS
  





! ******************************************************************************
!
! Purpose: Check whether geometry is correct for 1d and 2d runs.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_123d_checkgeometrywrapper(pRegion)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Parameters      
! ==============================================================================

    TYPE(t_region), POINTER :: pregion

! ==============================================================================
!   Locals      
! ==============================================================================

    TYPE(t_global), POINTER :: global      
    TYPE(t_grid), POINTER :: pgrid

! ******************************************************************************
!   Start
! ******************************************************************************

    global => pregion%global

    CALL registerfunction(global,'RFLU_123D_CheckGeometryWrapper',&
  'RFLU_ModDimensionality.F90')

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

    pgrid => pregion%grid

! ******************************************************************************
!   Check geometry for 2d runs
! ******************************************************************************

    SELECT CASE ( pregion%mixtInput%dimens ) 
      CASE ( 1 ) 
        CALL rflu_123d_checkgeometrykernel(pregion,ycoord)
        CALL rflu_123d_checkgeometrykernel(pregion,zcoord)      
      CASE ( 2 ) 
        CALL rflu_123d_checkgeometrykernel(pregion,zcoord)
      CASE ( 3 ) ! Defensive coding
      CASE default
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens

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

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_123d_checkgeometrywrapper








! ******************************************************************************
!
! Purpose: Check whether geometry is correct for 1d and 2d runs.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!   dir         Coordinate direction
!
! Output: None.
!
! Notes: 
!   1. dir-component of normals MUST be close to machine precision for all 
!      patches except virtual ones, for which they must be close to unity.
!
! ******************************************************************************

  SUBROUTINE rflu_123d_checkgeometrykernel(pRegion,dir)

    USE modtools, ONLY: floatequal

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Parameters      
! ==============================================================================

    INTEGER, INTENT(IN) :: dir
    TYPE(t_region), POINTER :: pregion

! ==============================================================================
!   Locals      
! ==============================================================================

    INTEGER :: ipatch
    REAL(RFREAL) :: ndmax,nxmax,ndmin,nxmin,ntol
    TYPE(t_grid), POINTER :: pgrid      
    TYPE(t_patch), POINTER :: ppatch            
    TYPE(t_global), POINTER :: global      

! ******************************************************************************
!   Start
! ******************************************************************************

    global => pregion%global

    CALL registerfunction(global,'RFLU_123D_CheckGeometryKernel',&
  'RFLU_ModDimensionality.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel >= verbose_high) THEN
      WRITE(stdout,'(A,1X,A)') solver_name,'Checking geometry...'
      WRITE(stdout,'(A,3X,A,1X,I2)') solver_name,'Component:',dir
    END IF ! global%myProcid

    pgrid => pregion%grid

    ntol =  1.0e-14_rfreal

! ******************************************************************************
!   Interior faces: dir-component of normal must be close to zero
! ******************************************************************************

    IF ( pgrid%nFacesTot > 0 ) THEN 
      ndmin = minval(pgrid%fn(dir,1:pgrid%nFacesTot))
      ndmax = maxval(pgrid%fn(dir,1:pgrid%nFacesTot))      

      IF ( global%myProcid == masterproc .AND. & 
           global%verbLevel >= verbose_high) THEN     
        WRITE(stdout,'(A,3X,A)') solver_name,'Extrema of face-normal vectors:'
        WRITE(stdout,'(A,5X,A,1X,E23.16)') solver_name,'Tolerance:',ntol
        WRITE(stdout,'(A,7X,A,1X,2(1X,E23.16))') solver_name,'Interior:', &
                                                 ndmin,ndmax
      END IF ! global%myProcid

      IF ( (abs(ndmin) > ntol) .OR. (abs(ndmax) > ntol) ) THEN 
        CALL errorstop(global,err_face_normal_invalid,__line__)
      END IF ! ABS(ndMin)
  
      IF ( (floatequal(abs(ndmin),0.0_rfreal,ntol) .EQV. .false.) .OR. & 
           (floatequal(abs(ndmax),0.0_rfreal,ntol) .EQV. .false.) ) THEN 
        CALL errorstop(global,err_face_normal_invalid,__line__)
      END IF ! FloatEqual     
    END IF ! pGrid%nFacesTot

! ******************************************************************************
!   Boundary faces: dir-component of normal must be close to zero for faces on
!   non-virtual patches, otherwise must be close to unity.
! ******************************************************************************

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch)

      IF ( ppatch%nBFacesTot > 0 ) THEN 
        ndmin = minval(ppatch%fn(dir,1:ppatch%nBFacesTot))
        ndmax = maxval(ppatch%fn(dir,1:ppatch%nBFacesTot))  

        nxmin = minval(ppatch%fn(xcoord,1:ppatch%nBFacesTot))
        nxmax = maxval(ppatch%fn(xcoord,1:ppatch%nBFacesTot))  

        IF ( global%myProcid == masterproc .AND. & 
             global%verbLevel >= verbose_high) THEN
          WRITE(stdout,'(A,7X,A,1X,I2,A,1X,2(1X,E23.16))') & 
                solver_name,'Patch',ipatch,':',ndmin,ndmax
        END IF ! global%myProcid

! ==============================================================================
!       Non-virtual patches: dir-component of normal must be close to zero 
! ==============================================================================

        IF ( ppatch%bcType /= bc_virtual ) THEN        
          IF ( (floatequal(abs(ndmin),0.0_rfreal,ntol) .EQV. .false.) .OR. & 
               (floatequal(abs(ndmax),0.0_rfreal,ntol) .EQV. .false.) ) THEN 
            CALL errorstop(global,err_face_normal_invalid,__line__)
          END IF ! FloatEqual 
  
! ==============================================================================
!       Virtual patches: x-component of normal must be close to zero
! ==============================================================================

        ELSE 
          IF ( (floatequal(abs(nxmin),0.0_rfreal,ntol) .EQV. .false.) .OR. & 
               (floatequal(abs(nxmax),0.0_rfreal,ntol) .EQV. .false.) ) THEN 
            CALL errorstop(global,err_face_normal_invalid,__line__)
          END IF ! FloatEqual     
        END IF ! pPatch%bcType            
      END IF ! pPatch%nBFacesTot
    END DO ! iPatch

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel >= verbose_high) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name,'Checking geometry done.'
    END IF ! global%verbLevel

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_123d_checkgeometrykernel






! ******************************************************************************
!
! Purpose: Check whether topology is correct.
!
! Description: None.
!
! Input:
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************

  SUBROUTINE rflu_123d_checktopology(pRegion)

    IMPLICIT NONE

! ******************************************************************************
!   Declarations and definitions
! ******************************************************************************

! ==============================================================================
!   Parameters      
! ==============================================================================

    TYPE(t_region), POINTER :: pregion

! ==============================================================================
!   Locals      
! ==============================================================================

    INTEGER :: icl,ifg,ifl,ifl2,ipatch,ipatchcntr
    TYPE(t_grid), POINTER :: pgrid      
    TYPE(t_patch), POINTER :: ppatch            
    TYPE(t_global), POINTER :: global      

! ******************************************************************************
!   Start
! ******************************************************************************

    global => pregion%global

    CALL registerfunction(global,'RFLU_123D_CheckTopology',&
  'RFLU_ModDimensionality.F90')

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel >= verbose_high) THEN
      WRITE(stdout,'(A,1X,A)') solver_name,'Checking topology...'
    END IF ! global%verbLevel

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

    pgrid => pregion%grid

! ******************************************************************************
!   Check dimensionality
! ******************************************************************************

! ==============================================================================
!   Check that have only appropriate cells     
! ==============================================================================
    
    SELECT CASE ( pregion%mixtInput%dimens ) 
      CASE ( 1 ) 
        IF ( pgrid%nTetsTot /= 0 .OR. &
             pgrid%nPrisTot /= 0 .OR. &
             pgrid%nPyrsTot /= 0 ) THEN
          CALL errorstop(global,err_dimens_invalid,__line__)
        END IF ! pGrid%nTetsTot    
      CASE ( 2 ) 
        IF ( pgrid%nTetsTot /= 0 .OR. pgrid%nPyrsTot /= 0 ) THEN
          CALL errorstop(global,err_dimens_invalid,__line__)
        END IF ! pGrid%nTetsTot
      CASE ( 3 ) 
      CASE default ! Defensive coding
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens
    
! ==============================================================================
!   Must have specified number of virtual patches    
! ==============================================================================

    ipatchcntr = 0 

    DO ipatch = 1,pgrid%nPatches
      ppatch => pregion%patches(ipatch) 

      IF ( ppatch%bcType == bc_virtual ) THEN 
        ipatchcntr = ipatchcntr + 1
      END IF ! pPatch%bcType
    END DO ! iPatch

    SELECT CASE ( pregion%mixtInput%dimens )
      CASE ( 1 )         
        IF ( ipatchcntr /= 4 ) THEN 
          CALL errorstop(global,err_num_bc_virtual,__line__)        
        END IF ! iPatchCntr    
      CASE ( 2 )         
        IF ( ipatchcntr /= 2 ) THEN 
          CALL errorstop(global,err_num_bc_virtual,__line__)        
        END IF ! iPatchCntr
      CASE ( 3 )
        IF ( ipatchcntr /= 0 ) THEN 
          CALL errorstop(global,err_num_bc_virtual,__line__)        
        END IF ! iPatchCntr       
      CASE default ! Defensive coding
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens    

! ==============================================================================
!   Each cell must have specified number of boundary faces. NOTE require 
!   cell-to-face connectivity array. For prisms, boundary faces must be 
!   triangular faces. 
! ==============================================================================

    SELECT CASE ( pregion%mixtInput%dimens )     

! ------------------------------------------------------------------------------
!     One dimension
! ------------------------------------------------------------------------------      

      CASE ( 1 ) 
      
! ----- Hexahedra --------------------------------------------------------------
      
        DO icl = 1,pgrid%nHexsTot
          ipatchcntr = 0 
        
          DO ifl = 1,6
            ipatch = pgrid%hex2f(1,ifl,icl)
            
            IF ( ipatch > 0 ) THEN
              ppatch => pregion%patches(ipatch)
             
              IF ( ppatch%bcType == bc_virtual ) THEN            
                ipatchcntr = ipatchcntr + 1
              END IF ! pPatch%bcType
            END IF ! iPatch
          END DO ! ifl
          
          IF ( ipatchcntr /= 4 ) THEN 
            CALL errorstop(global,err_dimens_invalid,__line__)
          END IF ! iPatchCntr
        END DO ! icl

! ------------------------------------------------------------------------------
!     Two dimensions
! ------------------------------------------------------------------------------      

      CASE ( 2 ) 
      
! ----- Hexahedra --------------------------------------------------------------
      
        DO icl = 1,pgrid%nHexsTot
          ipatchcntr = 0 
        
          DO ifl = 1,6
            ipatch = pgrid%hex2f(1,ifl,icl)
            
            IF ( ipatch > 0 ) THEN
              ppatch => pregion%patches(ipatch)
             
              IF ( ppatch%bcType == bc_virtual ) THEN            
                ipatchcntr = ipatchcntr + 1
              END IF ! pPatch%bcType
            END IF ! iPatch
          END DO ! ifl
          
          IF ( ipatchcntr /= 2 ) THEN 
            CALL errorstop(global,err_dimens_invalid,__line__)
          END IF ! iPatchCntr
        END DO ! icl

! ----- Prisms -----------------------------------------------------------------
        
        DO icl = 1,pgrid%nPrisTot
          ipatchcntr = 0 
        
          DO ifl = 1,5
            ipatch = pgrid%pri2f(1,ifl,icl)
            ifl2   = pgrid%pri2f(2,ifl,icl)
            
            IF ( ipatch > 0 ) THEN
              ppatch => pregion%patches(ipatch)
             
              IF ( ppatch%bcType == bc_virtual ) THEN            
                ipatchcntr = ipatchcntr + 1
                                
                IF ( ppatch%bf2v(4,ifl2) /= vert_none ) THEN 
                  CALL errorstop(global,err_dimens_invalid,__line__)
                END IF ! pGrid%f2v                
              END IF ! pPatch%bcType
            END IF ! iPatch            
          END DO ! ifl
          
          IF ( ipatchcntr /= 2 ) THEN 
            CALL errorstop(global,err_dimens_invalid,__line__)
          END IF ! iPatchCntr
        END DO ! icl       

! ------------------------------------------------------------------------------
!     Three dimensions and default
! ------------------------------------------------------------------------------      

      CASE ( 3 ) 
      CASE default ! Defensive coding
        CALL errorstop(global,err_reached_default,__line__)
    END SELECT ! pRegion%mixtInput%dimens

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

    IF ( global%myProcid == masterproc .AND. &
         global%verbLevel >= verbose_high) THEN 
      WRITE(stdout,'(A,1X,A)') solver_name,'Checking topology done.'
    END IF ! global%verbLevel

    CALL deregisterfunction(global)  

  END SUBROUTINE rflu_123d_checktopology




! ******************************************************************************
! End
! ******************************************************************************
  
END MODULE rflu_moddimensionality


! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ModDimensionality.F90,v $
! Revision 1.9  2008/12/06 08:44:21  mtcampbe
! Updated license.
!
! Revision 1.8  2008/11/19 22:17:32  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.7  2007/03/07 03:20:05  haselbac
! Added IFs on nFacesTot to avoid problems with single-cell grid
!
! Revision 1.6  2007/02/27 13:03:55  haselbac
! Enabled 1d computations
!
! Revision 1.5  2005/12/11 15:55:03  haselbac
! Bug fix: Added missing IF on myProcid
!
! Revision 1.4  2005/11/09 01:23:05  haselbac
! Increased tolerance again, did not increase version number
!
! Revision 1.3  2005/11/09 01:18:20  haselbac
! Changed tolerance after finding problems with Manojs cylinder comp
!
! Revision 1.2  2005/11/04 14:06:48  haselbac
! Renamed existing routine, added new routine to check geom
! 
! ******************************************************************************