RFLU_SetMoveGridOptions.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: Set options for moving grid computations.
!
! Description: The options for imposing boundary conditions on the grid motion
!   are determined by whether the patches are flat and whether they are aligned
!   with a coordinate direction. If not flat, no boundary condition is (can be)
!   applied without a geometrical description of the surface. If the patch is
!   flat, Dirichlet conditions are set if the patch normal is perfectly aligned
!   with a coordinate direction, otherwise von Neumann conditions are applied.
!
! Input: 
!   pRegion     Region pointer
!
! Output: None.
!
! Notes: 
!   1. At present, the only option which is determined in this routine is 
!      how the boundary conditions are set for the grid motion, namely either
!      von Neumann or Dirichlet conditions. The von Neumann condition is 
!      set for those patches the normal to which is not aligned with a 
!      coordinate direction. The Dirichlet condition is set for patches the 
!      normal to which is aligned with a coordinate direction.
!   2. This routine must be called after the boundary vertex normals are 
!      constructed, and after the boundary condition file is read (so that the
!      grid motion options for motion and smoothing are set).
!
! ******************************************************************************
!
! $Id: RFLU_SetMoveGridOptions.F90,v 1.15 2008/12/06 08:44:30 mtcampbe Exp $
!
! Copyright: (c) 2002-2005 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_setmovegridoptions(pRegion)

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

  USE modtools  
    
  IMPLICIT NONE

! **************************************I****************************************
! Definitions and declarations
! ******************************************************************************

! ==============================================================================
! Arguments
! ==============================================================================

  TYPE(t_region), POINTER :: pregion

! ==============================================================================
! Locals
! ==============================================================================
 
  LOGICAL :: outsiderflag,xflatflag,yflatflag,zflatflag
  CHARACTER(CHRLEN) :: movebctypestring,rcsidentstring  
  INTEGER :: errorflag,ibv,ipatch
  INTEGER, DIMENSION(2) :: maxnormloc
  REAL(RFREAL), PARAMETER :: eqtol = 1.0e-6_rfreal
  REAL(RFREAL) :: term,xnorm,xnormmax,xnormmin,ynorm,ynormmax,ynormmin, &  
                  znorm,znormmax,znormmin
  REAL(RFREAL), DIMENSION(:,:), ALLOCATABLE :: globalvals,localvals  
  REAL(RFREAL), DIMENSION(:,:,:), ALLOCATABLE :: bvnext
  TYPE(t_global), POINTER :: global
  TYPE(t_patch), POINTER :: ppatch

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

  rcsidentstring = '$RCSfile: RFLU_SetMoveGridOptions.F90,v $ $Revision: 1.15 $'

  global => pregion%global

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

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel > verbose_none ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name,'Setting grid motion options...'
    WRITE(stdout,'(A,3X,A,1X,E15.8)') solver_name,'Equality tolerance:',eqtol
    WRITE(stdout,'(A,3X,A,3X,A,2X,A)') solver_name,'Local','Global','Option'
  END IF ! global%verbLevel

! ******************************************************************************
! Allocate and initialize temporary memory
! ******************************************************************************

  ALLOCATE(bvnext(min_val:max_val,xcoord:zcoord,global%nPatches),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'bvnExt')
  END IF ! global%error

  DO ipatch = 1,global%nPatches
    bvnext(min_val,xcoord:zcoord,ipatch) =  huge(1.0_rfreal)
    bvnext(max_val,xcoord:zcoord,ipatch) = -huge(1.0_rfreal)              
  END DO ! iPatch

  ALLOCATE(localvals(xcoord:zcoord,global%nPatches),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'localVals')
  END IF ! global%error 
  
  ALLOCATE(globalvals(xcoord:zcoord,global%nPatches),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'globalVals')
  END IF ! global%error 

! ******************************************************************************
! Compute local extrema of patch vertex-normal vectors
! ******************************************************************************
 
  DO ipatch = 1,pregion%grid%nPatches
    ppatch => pregion%patches(ipatch)

    xnormmin =  huge(1.0_rfreal)
    xnormmax = -huge(1.0_rfreal)
    ynormmin =  huge(1.0_rfreal)
    ynormmax = -huge(1.0_rfreal)
    znormmin =  huge(1.0_rfreal)
    znormmax = -huge(1.0_rfreal)    

    DO ibv = 1,ppatch%nBVertTot
      xnorm = ppatch%bvn(xcoord,ibv) 
      ynorm = ppatch%bvn(ycoord,ibv)
      znorm = ppatch%bvn(zcoord,ibv)

      IF ( xnorm < xnormmin ) THEN 
        xnormmin = xnorm
      ELSE IF ( xnorm > xnormmax ) THEN 
        xnormmax = xnorm
      END IF ! xNorm

      IF ( ynorm < ynormmin ) THEN 
        ynormmin = ynorm
      ELSE IF ( ynorm > ynormmax ) THEN 
        ynormmax = ynorm
      END IF ! yNorm

      IF ( znorm < znormmin ) THEN 
        znormmin = znorm
      ELSE IF ( znorm > znormmax ) THEN 
        znormmax = znorm
      END IF ! zNorm
    END DO ! ibv

    bvnext(min_val,xcoord,ppatch%iPatchGlobal) = xnormmin
    bvnext(min_val,ycoord,ppatch%iPatchGlobal) = ynormmin
    bvnext(min_val,zcoord,ppatch%iPatchGlobal) = znormmin       

    bvnext(max_val,xcoord,ppatch%iPatchGlobal) = xnormmax
    bvnext(max_val,ycoord,ppatch%iPatchGlobal) = ynormmax
    bvnext(max_val,zcoord,ppatch%iPatchGlobal) = znormmax        
  END DO ! iPatch

! ******************************************************************************
! Compute global extrema of patch normal vectors
! ******************************************************************************

  DO ipatch = 1,global%nPatches
    localvals(xcoord,ipatch) = bvnext(min_val,xcoord,ipatch)
    localvals(ycoord,ipatch) = bvnext(min_val,ycoord,ipatch)
    localvals(zcoord,ipatch) = bvnext(min_val,zcoord,ipatch) 
  END DO ! iPatch  

  CALL mpi_reduce(localvals,globalvals,(zcoord-xcoord+1)*global%nPatches, & 
                  mpi_rfreal,mpi_min,masterproc,global%mpiComm,errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_mpi_trouble,__line__)
  END IF ! global%errorFlag
  
  DO ipatch = 1,global%nPatches
    bvnext(min_val,xcoord,ipatch) = globalvals(xcoord,ipatch)
    bvnext(min_val,ycoord,ipatch) = globalvals(ycoord,ipatch)
    bvnext(min_val,zcoord,ipatch) = globalvals(zcoord,ipatch)
        
    localvals(xcoord,ipatch) = bvnext(max_val,xcoord,ipatch)
    localvals(ycoord,ipatch) = bvnext(max_val,ycoord,ipatch)
    localvals(zcoord,ipatch) = bvnext(max_val,zcoord,ipatch)                      
  END DO ! iPatch  
   
  CALL mpi_reduce(localvals,globalvals,(zcoord-xcoord+1)*global%nPatches, & 
                  mpi_rfreal,mpi_min,masterproc,global%mpiComm,errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_mpi_trouble,__line__)
  END IF ! global%errorFlag
  
  DO ipatch = 1,global%nPatches
    bvnext(max_val,xcoord,ipatch) = globalvals(xcoord,ipatch)
    bvnext(max_val,ycoord,ipatch) = globalvals(ycoord,ipatch)
    bvnext(max_val,zcoord,ipatch) = globalvals(zcoord,ipatch)   
  END DO ! iPatch   

! ******************************************************************************
! Loop over patches to determine options on grid motion
! ******************************************************************************

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

    ppatch%moveBcType = movegrid_bctype_none  

! ==============================================================================
!   Determine whether non-moving patches are flat
! ==============================================================================

    IF ( (ppatch%movePatch .EQV. .false.) .AND. & 
         (ppatch%smoothGrid .EQV. .true.) ) THEN 

! ------------------------------------------------------------------------------
!     Get extrema of normal vector
! ------------------------------------------------------------------------------
    
      xnormmin = bvnext(min_val,xcoord,ppatch%iPatchGlobal)
      xnormmax = bvnext(max_val,xcoord,ppatch%iPatchGlobal)

      ynormmin = bvnext(min_val,ycoord,ppatch%iPatchGlobal)        
      ynormmax = bvnext(max_val,ycoord,ppatch%iPatchGlobal)

      znormmin = bvnext(min_val,zcoord,ppatch%iPatchGlobal)
      znormmax = bvnext(max_val,zcoord,ppatch%iPatchGlobal)
    
! ------------------------------------------------------------------------------
!     Determine whether patch is flat by checking whether extrema have same 
!     sign (if yes, and extrema are equal, have flat patch; if no, and extrema 
!     are close to zero, also have flat patch)
! ------------------------------------------------------------------------------
                
! --- x-direction --------------------------------------------------------------    
        
      IF ( (sign(1.0_rfreal,xnormmin) == sign(1.0_rfreal,xnormmax)) ) THEN 
        IF ( floatequal(xnormmin,xnormmax,eqtol) .EQV. .true. ) THEN 
          xflatflag = .true.
        ELSE 
          xflatflag = .false.
        END IF ! FloatEqual
      ELSE
         IF ( (floatequal(abs(xnormmin),epsilon(1.0_rfreal),eqtol)) .AND. & 
              (floatequal(abs(xnormmax),epsilon(1.0_rfreal),eqtol)) ) THEN
          xflatflag = .true.      
        ELSE
          xflatflag = .false.
        END IF ! FloatEqual
      END IF ! FloatEqual

! --- y-direction --------------------------------------------------------------    

      IF ( (sign(1.0_rfreal,ynormmin) == sign(1.0_rfreal,ynormmax)) ) THEN 
        IF ( floatequal(ynormmin,ynormmax,eqtol) .EQV. .true. ) THEN 
          yflatflag = .true.
        ELSE 
          yflatflag = .false.
        END IF ! FloatEqual
      ELSE
         IF ( (floatequal(abs(ynormmin),epsilon(1.0_rfreal),eqtol)) .AND. & 
              (floatequal(abs(ynormmax),epsilon(1.0_rfreal),eqtol)) ) THEN
          yflatflag = .true.      
        ELSE
          yflatflag = .false.
        END IF ! FloatEqual
      END IF ! FloatEqual

! --- z-direction --------------------------------------------------------------    

      IF ( (sign(1.0_rfreal,znormmin) == sign(1.0_rfreal,znormmax)) ) THEN     
        IF ( floatequal(znormmin,znormmax,eqtol) .EQV. .true. ) THEN 
          zflatflag = .true.
        ELSE 
          zflatflag = .false.
        END IF ! FloatEqual
      ELSE    
         IF ( (floatequal(abs(znormmin),epsilon(1.0_rfreal),eqtol)) .AND. & 
              (floatequal(abs(znormmax),epsilon(1.0_rfreal),eqtol)) ) THEN
          zflatflag = .true.      
        ELSE
          zflatflag = .false.
        END IF ! FloatEqual
      END IF ! FloatEqual
    
! ==============================================================================
!     Determine whether patch is flat and aligned with coordinate directions. 
!     If flat and aligned with coordinate directions, use Dirichlet conditions,
!     otherwise use Neumann conditions (but only if do not have any outsider
!     actual vertices).
! ==============================================================================
                  
      IF ( (xflatflag .EQV. .true.) .AND. & 
           (yflatflag .EQV. .true.) .AND. & 
           (zflatflag .EQV. .true.) ) THEN 
        term = max(abs(xnormmin),abs(xnormmax), & 
                   abs(ynormmin),abs(ynormmax), &
                   abs(znormmin),abs(znormmax))             
                   
! ------------------------------------------------------------------------------
!       Patch is flat and aligned with coordinate axes
! ------------------------------------------------------------------------------
        
        IF ( floatequal(term,1.0_rfreal) .EQV. .true. ) THEN  
          IF ( term == abs(xnormmin) .OR. term == abs(xnormmax) ) THEN
            ppatch%moveBcType = xcoord       
          ELSE IF ( term == abs(ynormmin) .OR. term == abs(ynormmax) ) THEN
            ppatch%moveBcType = ycoord
          ELSE IF ( term == abs(znormmin) .OR. term == abs(znormmax) ) THEN
            ppatch%moveBcType = zcoord  
          END IF ! term                                  

! ------------------------------------------------------------------------------
!       Patch is flat but not aligned with coordinate axes, so set boundary 
!       condition to von Neumann
! ------------------------------------------------------------------------------        
        
        ELSE
          ppatch%moveBcType = movegrid_bctype_neumann
        END IF ! xFlagFlag
      END IF ! xFlatFlag

! ==============================================================================
!     Check that patches have proper boundary condition. Patches with active 
!     smoothing must be flat, whether or not aligned with coordinate directions.
! ==============================================================================

      IF ( ppatch%moveBcType == movegrid_bctype_none ) THEN 
        ppatch%smoothGrid = .false.
         
        global%warnCounter = global%warnCounter + 1  
         
        IF ( global%myProcid == masterproc .AND. & 
             global%verbLevel > verbose_none ) THEN            
           WRITE(stdout,'(A,3X,A,I3,A)') solver_name, &
                 '*** WARNING *** Invalid smoothing input for patch ',ipatch, &
                 '. Overriding user input for grid smoothing.'
        END IF ! global%myProcid                      
      END IF ! pPatch%smoothGrid
    END IF ! pPatch%movePatch

! ==============================================================================
!   Write out info
! ==============================================================================
  
    IF ( global%myProcid == masterproc .AND. & 
         global%verbLevel > verbose_none ) THEN
      IF ( ppatch%moveBcType == movegrid_bctype_none ) THEN 
        movebctypestring = 'None'   
      ELSE IF ( ppatch%moveBcType == movegrid_bctype_neumann ) THEN 
        movebctypestring = 'von Neumann (homogeneous)'         
      ELSE IF ( ppatch%moveBcType == movegrid_bctype_dirichx ) THEN 
        movebctypestring = 'Dirichlet (homogeneous in x)'     
      ELSE IF ( ppatch%moveBcType == movegrid_bctype_dirichy ) THEN 
        movebctypestring = 'Dirichlet (homogeneous in y)' 
      ELSE IF ( ppatch%moveBcType == movegrid_bctype_dirichz ) THEN 
        movebctypestring = 'Dirichlet (homogeneous in z)' 
      ELSE 
        CALL errorstop(global,err_reached_default,__line__)
      END IF ! pPatch%moveBcType          
           
      WRITE(stdout,'(A,2X,I4,5X,I4,4X,A)') solver_name,ipatch, & 
                                           ppatch%iPatchGlobal, & 
                                           trim(movebctypestring)
    END IF ! global%myProcid
  END DO ! iPatch

! ******************************************************************************
! Check settings - catch error before get to grid motion (caught there also)
! ******************************************************************************

  DO ipatch = 1,pregion%grid%nPatches
    ppatch => pregion%patches(ipatch)
  
    IF ( (ppatch%smoothGrid .EQV. .true.) .AND. & 
         (ppatch%moveBcType == movegrid_bctype_none) ) THEN 
      CALL errorstop(global,err_movepatch_bc_invalid,__line__)
    END IF ! pPatch
  END DO ! iPatch

! ******************************************************************************
! Deallocate temporary memory
! ******************************************************************************

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

  DEALLOCATE(localvals,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'localVals')
  END IF ! global%error 
  
  DEALLOCATE(globalvals,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'globalVals')
  END IF ! global%error 

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

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel > verbose_none ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name,'Setting grid motion options done.'
  END IF ! global%verbLevel

  CALL deregisterfunction(global)

END SUBROUTINE rflu_setmovegridoptions

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_SetMoveGridOptions.F90,v $
! Revision 1.15  2008/12/06 08:44:30  mtcampbe
! Updated license.
!
! Revision 1.14  2008/11/19 22:17:43  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.13  2006/04/07 15:19:22  haselbac
! Removed tabs
!
! Revision 1.12  2005/04/15 15:07:25  haselbac
! Converted to MPI
!
! Revision 1.11  2004/10/19 19:29:24  haselbac
! Removed logic related to insider/outsider vertices
!
! Revision 1.10  2003/07/22 02:12:00  haselbac
! Added global%warnCounter
!
! Revision 1.9  2003/05/01 14:13:51  haselbac
! Added logic to deal with non-moving and non-smoothed patches
!
! Revision 1.8  2003/03/27 21:50:21  haselbac
! Fixed bug in determining outsiderFlag
!
! Revision 1.7  2003/03/15 18:57:22  haselbac
! Completed || of gm, added check for outsider vertices
!
! Revision 1.6  2003/02/20 20:15:23  haselbac
! Bug fix: nVertTot changed to nVert (from different working version)
!
! Revision 1.5  2003/02/20 19:49:20  haselbac
! Corrected bug in check
!
! Revision 1.4  2003/02/06 19:32:16  haselbac
! Bug fix (EPS instd of PREC), use tolerance, check for error
!
! Revision 1.3  2003/01/28 14:49:19  haselbac
! Only set bc for non-moving patches
!
! Revision 1.2  2002/11/27 20:27:17  haselbac
! Changed test for planarity and added output
!
! Revision 1.1  2002/11/26 15:28:56  haselbac
! Initial revision
!
! ******************************************************************************