RFLU_MoveGridXyz.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: Move grid by smoothing coordinates.
!
! Description: Simple pseudo-Laplacian smoothing of grid given boundary 
!   distortion. 
!
! Input: 
!   regions     Region data
!   context     Context in which grid smoothing is to be done 
!
! Output: None.
!
! Notes: 
!   1. The input variable context allows the routine to distinguish between 
!      two kinds of grid smoothing operations: The first is an operation on 
!      a grid with the goal of smoothing the interior grid with no boundary
!      motion. This could be done to smooth an original grid. The second is 
!      the smoothing of the grid with boundary motion. Since the operations 
!      to be carried out are slighly different, it is convenient to use the
!      context variable to distinguish between them. The context variable can
!      take the two (integer parameter) values MOVEGRID_CONTEXT_ONLYSMOOTH 
!      and MOVEGRID_CONTEXT_MOVESMOOTH.
!   2. The two parameters influencing the grid motion, i.e., the number of
!      iterations and the smoothing coefficient, may need some tuning as more
!      complicated geometries are attempted...
!   3. The integrity of patches with no imposed motion is ensured by one of 
!      two methods. Both may be cast in terms of a boundary condition on the 
!      residual of the Laplacian operator. The first is a von Neumann bc on 
!      the residual and can be used for flat boundaries which are not aligned
!      the normal to which is not aligned with a coordinate direction. The 
!      second is a Dirichlet boundary condition which can be used for flat 
!      boundaries the normal to which is aligned with a coordinate direction.
!      Of course, the first method subsumes the second, but the latter may be
!      preferable because the first can introduce errors on the order of the 
!      machine precision into the constant coordinate component. Which of the 
!      two methods is more suitable is determined by a simple algorithm in
!      the routine RFLU_SetMoveGridOptions.F90. 
!
! ******************************************************************************
!
! $Id: RFLU_MoveGridXyz.F90,v 1.9 2008/12/06 08:44:30 mtcampbe Exp $
!
! Copyright: (c) 2002-2005 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_movegridxyz(regions,context)

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

  IMPLICIT NONE

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

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

  INTEGER, INTENT(IN) :: context
  TYPE(t_region), DIMENSION(:), POINTER :: regions

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

  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER :: errorflag,ic,ie,ipatch,ireg,it,iv,ivg,niter,patchcheckcounter, & 
             v1,v2,v3,v4
  REAL(RFREAL) :: nx,ny,nz,sfact,term,volmaxglob,volmaxloc,volminglob, &
                  volminloc
  REAL(RFREAL) :: dxyz(xcoord:zcoord)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: prhs,pxyz,pxyzold
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid,pgridold
  TYPE(t_patch), POINTER :: ppatch
  TYPE(t_region), POINTER :: pregion

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

  rcsidentstring = '$RCSfile: RFLU_MoveGridXyz.F90,v $ $Revision: 1.9 $'

  global => regions(1)%global

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

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel /= verbose_low ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name,'Moving grid based on coordinates...'
  END IF ! global%myProcid
  
! ******************************************************************************
! Set variables
! ******************************************************************************
    
  niter = regions(1)%mixtInput%moveGridNIter
  sfact = regions(1)%mixtInput%moveGridSFact 

! TEMPORARY Disable smoothing for MPI version
  niter = 0 
! END TEMPORARY
  
! ******************************************************************************
! Copy grid data to old grid
! ******************************************************************************

  DO ireg = 1,global%nRegionsLocal 
    pgrid    => regions(ireg)%grid
    pgridold => regions(ireg)%gridOld    
    
    IF ( context == movegrid_context_movesmooth ) THEN 
      DO ic = 1,pgrid%nCellsTot ! Explicit copy to avoid ASCI White problem
        pgridold%vol(ic) = pgrid%vol(ic)
      END DO ! ic  
    END IF ! context   
         
    DO iv = 1,pgrid%nVertTot ! Explicit copy to avoid ASCI White problem    
      pgridold%xyz(xcoord,iv) = pgrid%xyz(xcoord,iv)
      pgridold%xyz(ycoord,iv) = pgrid%xyz(ycoord,iv)
      pgridold%xyz(zcoord,iv) = pgrid%xyz(zcoord,iv)            
    END DO ! iv 
  END DO ! iReg
  
! ******************************************************************************
! Allocate vertex degree array
! ******************************************************************************  
  
  DO ireg = 1,global%nRegionsLocal
    pgrid => regions(ireg)%grid
       
    ALLOCATE(pgrid%degr(pgrid%nVertTot),stat=errorflag)
    global%error = errorflag
    IF ( global%error /= err_none ) THEN 
      CALL errorstop(global,err_allocate,__line__,'pGrid%degr')
    END IF ! global%error  
        
    IF ( pgrid%nTetsTot == pgrid%nCellsTot ) THEN ! Only for tetrahedra
      ALLOCATE(pgrid%volMin(pgrid%nVertTot),stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN
        CALL errorstop(global,err_allocate,__line__,'pGrid%volMin')
      END IF ! global%error 
    END IF ! pGrid%nTetsTot    
  END DO ! iReg  
    
! ******************************************************************************  
! Compute modified edge degree. NOTE the modified vertex degree depends on the
! geometric properties of the grid; in this particular version on the volume. 
! Because the volume does not get updated during the iterations to smooth the
! grid, the modified edge degree does not have to be recomputed. 
! ******************************************************************************  
  
! ==============================================================================  
! Original vertex degree. NOTE must take into account virtual edges.
! ==============================================================================  

  DO ireg = 1,global%nRegionsLocal
    pgrid => regions(ireg)%grid

    DO iv = 1,pgrid%nVertTot ! Explicit loop to avoid ASCI White problem 
      pgrid%degr(iv) = 0
    END DO ! iv 

    DO ie = 1,pgrid%nEdgesTot
      v1 = pgrid%e2v(1,ie)
      v2 = pgrid%e2v(2,ie)

      pgrid%degr(v1) = pgrid%degr(v1) + 1
      pgrid%degr(v2) = pgrid%degr(v2) + 1      
    END DO ! ie
  END DO ! iReg     

! ==============================================================================  
! Modified vertex degree (for tetrahedral grids only). NOTE the modification 
! is usually beneficial if the grid contains very small and very large cells 
! and the small cells are of poor quality, which can lead to negative volumes.
! ==============================================================================  
        
  DO ireg = 1,global%nRegionsLocal 
    pgrid => regions(ireg)%grid

    IF ( pgrid%nTetsTot == pgrid%nCellsTot ) THEN 
      volminloc = minval(pgrid%vol(1:pgrid%nTets))
      volmaxloc = maxval(pgrid%vol(1:pgrid%nTets))

! TO DO 
! Need to use reductions to find volume extrema
      volminglob = volminloc
      volmaxglob = volmaxloc  
! END TO DO 
   
      DO iv = 1,pgrid%nVertTot ! Explicit loop to avoid ASCI White problem 
        pgrid%volMin(iv) = huge(1.0_rfreal)
      END DO ! iv

      DO ic = 1,pgrid%nTetsTot
        v1 = pgrid%tet2v(1,ic)
        v2 = pgrid%tet2v(2,ic)
        v3 = pgrid%tet2v(3,ic)
        v4 = pgrid%tet2v(4,ic) 

        pgrid%volMin(v1) = min(pgrid%volMin(v1),pgrid%vol(ic))
        pgrid%volMin(v2) = min(pgrid%volMin(v2),pgrid%vol(ic))
        pgrid%volMin(v3) = min(pgrid%volMin(v3),pgrid%vol(ic))            
        pgrid%volMin(v4) = min(pgrid%volMin(v4),pgrid%vol(ic))
      END DO ! ic

      DO iv = 1,pgrid%nVert
        pgrid%degr(iv) = pgrid%degr(iv) & 
                       + (volmaxglob - volminglob)/pgrid%volMin(iv)
      END DO ! iv        
    END IF ! pGrid%nTetsTot
  END DO ! iReg       
    
! ******************************************************************************
! Initialize coordinates to old coordinates plus patch deformation
! ******************************************************************************
  
  DO ireg = 1,global%nRegionsLocal 
    pgrid    => regions(ireg)%grid
    pgridold => regions(ireg)%gridOld     
  
    pxyz    => pgrid%xyz
    pxyzold => pgridold%xyz    
    
    IF ( context == movegrid_context_movesmooth ) THEN 
      DO ipatch = 1,pgrid%nPatches
        ppatch => regions(ireg)%patches(ipatch)

        IF ( ppatch%movePatchDir /= movepatch_dir_none ) THEN       
          DO iv = 1,ppatch%nBVert
            ivg = ppatch%bv(iv)            

            pxyz(xcoord,ivg) = pxyzold(xcoord,ivg) + ppatch%dXyz(xcoord,iv)
            pxyz(ycoord,ivg) = pxyzold(ycoord,ivg) + ppatch%dXyz(ycoord,iv)
            pxyz(zcoord,ivg) = pxyzold(zcoord,ivg) + ppatch%dXyz(zcoord,iv)
          END DO ! iv
        END IF ! pPatch%movePatchDir
      END DO ! iPatch
    ELSE IF ( context == movegrid_context_onlysmooth ) THEN 
      DO ipatch = 1,pgrid%nPatches
        ppatch => regions(ireg)%patches(ipatch)

        IF ( ppatch%movePatchDir /= movepatch_dir_none ) THEN       
          DO iv = 1,ppatch%nBVert
            ivg = ppatch%bv(iv)

            pxyz(xcoord,ivg) = pxyzold(xcoord,ivg)
            pxyz(ycoord,ivg) = pxyzold(ycoord,ivg)
            pxyz(zcoord,ivg) = pxyzold(zcoord,ivg)
          END DO ! iv
        END IF ! pPatch%movePatchDir
      END DO ! iPatch    
    END IF ! context
  END DO ! iReg      
     
! ******************************************************************************
! Apply pseudo-Laplacian smoothing to coordinates of all vertices
! ******************************************************************************

  DO it = 1,niter
             
! ==============================================================================
!   Compute right-hand side with boundary restrictions already imposed
! ==============================================================================  
   
    DO ireg = 1,global%nRegionsLocal  
      pgrid => regions(ireg)%grid
      pxyz  => pgrid%xyz
      prhs  => pgrid%rhs     
      
! ------------------------------------------------------------------------------
!     Initialize right-hand side
! ------------------------------------------------------------------------------
      
      DO iv = 1,pgrid%nVertTot ! Explicit loop to avoid ASCI White problem 
        prhs(xcoord,iv) = 0.0_rfreal
        prhs(ycoord,iv) = 0.0_rfreal
        prhs(zcoord,iv) = 0.0_rfreal                
      END DO ! iv        
      
! ------------------------------------------------------------------------------
!     Compute right-hand side. NOTE loop only over actual edges and take into
!     accont the region degree of each edge to deal with multiple contributions
!     of edges on region boundaries.
! ------------------------------------------------------------------------------
  
      DO ie = 1,pgrid%nEdges 
        v1 = pgrid%e2v(1,ie)
        v2 = pgrid%e2v(2,ie)

        term = 1.0_rfreal/(REAL(pGrid%e2rDegr(ie),kind=rfreal))

        dxyz(xcoord) = term*(pxyz(xcoord,v2) - pxyz(xcoord,v1))
        dxyz(ycoord) = term*(pxyz(ycoord,v2) - pxyz(ycoord,v1))
        dxyz(zcoord) = term*(pxyz(zcoord,v2) - pxyz(zcoord,v1)) 

        prhs(xcoord,v1) = prhs(xcoord,v1) + dxyz(xcoord)
        prhs(ycoord,v1) = prhs(ycoord,v1) + dxyz(ycoord)
        prhs(zcoord,v1) = prhs(zcoord,v1) + dxyz(zcoord)     

        prhs(xcoord,v2) = prhs(xcoord,v2) - dxyz(xcoord)
        prhs(ycoord,v2) = prhs(ycoord,v2) - dxyz(ycoord)
        prhs(zcoord,v2) = prhs(zcoord,v2) - dxyz(zcoord)         
      END DO ! ie

! ------------------------------------------------------------------------------
!     Apply boundary deformation restrictions
! ------------------------------------------------------------------------------  
  
      patchcheckcounter = 0
  
      IF ( context == movegrid_context_movesmooth ) THEN 
      
! ----- Enforce imposed deformation --------------------------------------------
      
        DO ipatch = 1,pgrid%nPatches
          ppatch => regions(ireg)%patches(ipatch)

          IF ( ppatch%movePatchDir /= movepatch_dir_none ) THEN                  
            patchcheckcounter = patchcheckcounter + 1   
                   
            DO iv = 1,ppatch%nBVert
              ivg = ppatch%bv(iv)

              prhs(xcoord,ivg) = 0.0_rfreal
              prhs(ycoord,ivg) = 0.0_rfreal
              prhs(zcoord,ivg) = 0.0_rfreal
            END DO ! iv
          END IF ! pPatch%movePatch
        END DO ! iPatch                   
        
! ----- Ensure boundary integrity (no imposed motion or smoothing) -------------   
        
        DO ipatch = 1,pgrid%nPatches
          ppatch => regions(ireg)%patches(ipatch)

          IF ( ppatch%movePatchDir == movepatch_dir_none ) THEN                 
            patchcheckcounter = patchcheckcounter + 1                  
                 
            SELECT CASE ( ppatch%moveBcType ) 
            
! ----------- Normal displacement set to zero with Neumann bc              
            
              CASE ( movegrid_bctype_neumann )
                DO iv = 1,ppatch%nBVert
                  ivg = ppatch%bv(iv)

                  nx = ppatch%bvn(xcoord,iv)
                  ny = ppatch%bvn(ycoord,iv)
                  nz = ppatch%bvn(zcoord,iv)                    

                  term = prhs(xcoord,ivg)*nx & 
                       + prhs(ycoord,ivg)*ny & 
                       + prhs(zcoord,ivg)*nz

                  prhs(xcoord,ivg) = prhs(xcoord,ivg) - term*nx
                  prhs(ycoord,ivg) = prhs(ycoord,ivg) - term*ny
                  prhs(zcoord,ivg) = prhs(zcoord,ivg) - term*nz            
                END DO ! iv
                
! ----------- Normal displacement set to zero with Dirichlet bc                 
                
              CASE ( movegrid_bctype_dirichx, & 
                     movegrid_bctype_dirichy, & 
                     movegrid_bctype_dirichz )
                DO iv = 1,ppatch%nBVert
                  ivg = ppatch%bv(iv)
                  
                  prhs(ppatch%moveBcType,ivg) = 0.0_rfreal               
                END DO ! iv

! ----------- No bc (must only occur with no movement and smoothing)

              CASE ( movegrid_bctype_none ) 
                IF ( ppatch%smoothGrid .EQV. .false. ) THEN 
                  DO iv = 1,ppatch%nBVert
                    ivg = ppatch%bv(iv)

                    prhs(xcoord,ivg) = 0.0_rfreal
                    prhs(ycoord,ivg) = 0.0_rfreal
                    prhs(zcoord,ivg) = 0.0_rfreal
                  END DO ! iv                  
                ELSE ! Defensive programming
                  CALL errorstop(global,err_reached_default,__line__)
                END IF ! 

              CASE default ! Defensive programming
                CALL errorstop(global,err_reached_default,__line__)
            END SELECT ! pPatch%moveBcType 
                      
          END IF ! pPatch
        END DO ! iPatch
        
      ELSE IF ( context == movegrid_context_onlysmooth ) THEN 
        DO ipatch = 1,pgrid%nPatches
          ppatch => regions(ireg)%patches(ipatch)
       
          patchcheckcounter = patchcheckcounter + 1 
       
          DO iv = 1,ppatch%nBVert
            ivg = ppatch%bv(iv)

            prhs(xcoord,ivg) = 0.0_rfreal
            prhs(ycoord,ivg) = 0.0_rfreal
            prhs(zcoord,ivg) = 0.0_rfreal
          END DO ! iv
        END DO ! iPatch      
      END IF ! context
    END DO ! iReg

! ------------------------------------------------------------------------------
!   Check that all patches have mesh motion boundary conditions set. This 
!   should never happen - if it does, there is a bug in the above code 
!   segment (defensive coding)
! ------------------------------------------------------------------------------  
            
    IF ( patchcheckcounter /= pgrid%nPatches ) THEN 
      CALL errorstop(global,err_movepatch_bc_notset,__line__)
    END IF ! patchCheckCounter   

! TO DO     
! Need to use reduction to add contributions at interface vertices 
! END TO DO 
              
! ==============================================================================
!   Update coordinates of actual vertices
! ==============================================================================

    DO ireg = 1,global%nRegionsLocal
      pgrid => regions(ireg)%grid
      pxyz  => pgrid%xyz
      prhs  => pgrid%rhs    
    
      DO iv = 1,pgrid%nVert
        term = sfact/REAL(pGrid%degr(iv),RFREAL)     
      
        pxyz(xcoord,iv) = pxyz(xcoord,iv) + term*prhs(xcoord,iv)
        pxyz(ycoord,iv) = pxyz(ycoord,iv) + term*prhs(ycoord,iv)
        pxyz(zcoord,iv) = pxyz(zcoord,iv) + term*prhs(zcoord,iv)
      END DO ! iv
    END DO ! iReg

! TO DO 
! Update virtual vertices
! END TO DO 
  END DO ! it

! ******************************************************************************
! Deallocate residual and degree arrays
! ******************************************************************************  
  
  DO ireg = 1,global%nRegionsLocal 
    pgrid => regions(ireg)%grid
        
    DEALLOCATE(pgrid%degr,stat=errorflag)
    global%error = errorflag
    IF ( global%error /= err_none ) THEN 
      CALL errorstop(global,err_deallocate,__line__,'pGrid%degr')
    END IF ! global%error 
    
    IF ( pgrid%nTetsTot == pgrid%nCellsTot ) THEN ! Only for tetrahedra
      DEALLOCATE(pgrid%volMin,stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN
        CALL errorstop(global,err_deallocate,__line__,'pGrid%volMin')
      END IF ! global%error 
    END IF ! pGrid%nTetsTot             
  END DO ! iReg  

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

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel /= verbose_low ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name, &
                             'Moving grid based on coordinates done.'
  END IF ! global%myProcid

  CALL deregisterfunction(global)

END SUBROUTINE rflu_movegridxyz

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_MoveGridXyz.F90,v $
! Revision 1.9  2008/12/06 08:44:30  mtcampbe
! Updated license.
!
! Revision 1.8  2008/11/19 22:17:42  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.7  2006/04/07 15:19:22  haselbac
! Removed tabs
!
! Revision 1.6  2005/06/09 20:22:58  haselbac
! Replaced movePatch by movePatchDir
!
! Revision 1.5  2005/04/15 15:07:21  haselbac
! Converted to MPI, grid motion can only be used for serial runs
!
! Revision 1.4  2003/08/20 02:09:58  haselbac
! Changed verbosity conditions to reduce solver output in GENx runs
!
! Revision 1.3  2003/07/09 22:36:53  haselbac
! Added IF statement to avoid probs with single-proc charm jobs
!
! Revision 1.2  2003/05/01 14:12:03  haselbac
! Added logic to deal with non-moving and non-smoothed patches
!
! Revision 1.1  2003/03/31 16:05:39  haselbac
! Initial revision
!
! Revision 1.8  2003/03/27 19:31:36  haselbac
! Fixed bug in initializing pGrid%degr
!
! Revision 1.7  2003/03/15 18:51:03  haselbac
! Completed || of grid motion
!
! Revision 1.6  2003/02/20 19:43:51  haselbac
! Some simplification, added additional check
!
! Revision 1.5  2003/01/28 14:43:19  haselbac
! Now get user input, enforce no smoothing if necessary, small area modification activated
!
! Revision 1.4  2002/11/26 15:28:35  haselbac
! Added Dirichlet bc for grid motion
!
! Revision 1.3  2002/11/19 23:46:41  haselbac
! Cleaned up boundary integrity enforcement of grid motion
!
! Revision 1.2  2002/11/08 21:33:11  haselbac
! Several modifications, mainly cosmetic apart from degree-modification
!
! Revision 1.1  2002/10/27 19:15:20  haselbac
! Initial revision
!
! Revision 1.1  2002/10/16 21:15:11  haselbac
! Initial revision
!
! ******************************************************************************