RFLU_MoveGridDisp.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 displacements.
!
! Description: Simple pseudo-Laplacian smoothing of grid given boundary 
!   distortion. 
!
! Input: 
!   regions     Region data
!
! Output: None.
!
! Notes: 
!   1. 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...
!   2. 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_MoveGridDisp.F90,v 1.11 2008/12/06 08:44:30 mtcampbe Exp $
!
! Copyright: (c) 2003-2005 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_movegriddisp(regions)

  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
! =============================================================================

  TYPE(t_region), DIMENSION(:), POINTER :: regions

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

  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER :: errorflag,ic,ie,ipatch,ireg,it,iv,ivg,niter,patchcheckcounter, & 
             v1,v2
  REAL(RFREAL) :: nx,ny,nz,sfact,term,x1,x2,y1,y2,z1,z2
  REAL(RFREAL) :: ddisp(xcoord:zcoord)
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pdisp,prhs,pxyz
  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_MoveGridDisp.F90,v $ $Revision: 1.11 $'

  global => regions(1)%global

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

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel /= verbose_low ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name, &
                             'Moving grid based on displacements...'
  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. Must be done even if have zero smoothing 
! iterations. 
! *****************************************************************************

  DO ireg = 1,global%nRegionsLocal 
    pgrid    => regions(ireg)%grid
    pgridold => regions(ireg)%gridOld    
    
    DO ic = 1,pgrid%nCellsTot ! Explicit copy to avoid ASCI White problem
      pgridold%vol(ic) = pgrid%vol(ic)
    END DO ! ic  
         
    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 temporary memory for and compute edge weights 
! *****************************************************************************  
    
  IF ( niter > 0 ) THEN 
    DO ireg = 1,global%nRegionsLocal
      pgrid => regions(ireg)%grid
       
      ALLOCATE(pgrid%gmEdgeWght(pgrid%nEdgesTot),stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_allocate,__line__,'pGrid%gmEdgeWght')
      END IF ! global%error     )   
       
      ALLOCATE(pgrid%gmVertWght(pgrid%nVertTot),stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_allocate,__line__,'pGrid%gmVertWght')
      END IF ! global%error     
    END DO ! iReg  
    
    DO ireg = 1,global%nRegionsLocal
      pgrid => regions(ireg)%grid
  
      DO iv = 1,pgrid%nVertTot
        pgrid%gmVertWght(iv) = 0.0_rfreal
      END DO ! iv

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

        x1 = pgrid%xyz(xcoord,v1)
        y1 = pgrid%xyz(ycoord,v1)
        z1 = pgrid%xyz(zcoord,v1)
      
        x2 = pgrid%xyz(xcoord,v2)
        y2 = pgrid%xyz(ycoord,v2)
        z2 = pgrid%xyz(zcoord,v2)                  

        term = 1.0_rfreal/sqrt((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)
  
        pgrid%gmEdgeWght(ie) = term
      
        pgrid%gmVertWght(v1) = pgrid%gmVertWght(v1) + term
        pgrid%gmVertWght(v2) = pgrid%gmVertWght(v2) + term    
      END DO ! ie    
    END DO ! iReg     
  END IF ! nIter
    
! ******************************************************************************
! Initialize displacements to zero and impose patch displacements. NOTE need to
! impose patch displacements so that get boundary deformation even if have zero
! smoothing iterations. 
! ******************************************************************************    
  
  DO ireg = 1,global%nRegionsLocal 
    pgrid => regions(ireg)%grid
    pdisp => pgrid%disp  
    
    DO iv = 1,pgrid%nVertTot
      pdisp(xcoord,iv) = 0.0_rfreal 
      pdisp(ycoord,iv) = 0.0_rfreal       
      pdisp(zcoord,iv) = 0.0_rfreal          
    END DO ! iv 
    
    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)            
 
          pdisp(xcoord,ivg) = ppatch%dXyz(xcoord,iv)
          pdisp(ycoord,ivg) = ppatch%dXyz(ycoord,iv)
          pdisp(zcoord,ivg) = ppatch%dXyz(zcoord,iv)
        END DO ! iv
      END IF ! pPatch%movePatchDir
    END DO ! iPatch         
  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
      pdisp => pgrid%disp 
      prhs  => pgrid%rhs           
      
      DO iv = 1,pgrid%nVert
        prhs(xcoord,iv) = 0.0_rfreal  
        prhs(ycoord,iv) = 0.0_rfreal              
        prhs(zcoord,iv) = 0.0_rfreal  
      END DO ! iv    
      
! ------------------------------------------------------------------------------
!     Enforce patch deformation. NOTE that strictly speaking this is not
!     needed during first smoothing iteration because initialized outside 
!     iteration loop. 
! ------------------------------------------------------------------------------      
            
      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)            

            pdisp(xcoord,ivg) = ppatch%dXyz(xcoord,iv)
            pdisp(ycoord,ivg) = ppatch%dXyz(ycoord,iv)
            pdisp(zcoord,ivg) = ppatch%dXyz(zcoord,iv)
          END DO ! iv
        END IF ! pPatch%movePatchDir
      END DO ! iPatch      
      
! ------------------------------------------------------------------------------
!     Compute right-hand side. NOTE loop only over actual edges and take into
!     account 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 = sfact*pgrid%gmEdgeWght(ie)/(REAL(pGrid%e2rDegr(ie),kind=rfreal))

        ddisp(xcoord) = term*(pdisp(xcoord,v2) - pdisp(xcoord,v1))
        ddisp(ycoord) = term*(pdisp(ycoord,v2) - pdisp(ycoord,v1))
        ddisp(zcoord) = term*(pdisp(zcoord,v2) - pdisp(zcoord,v1)) 

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

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

! ------------------------------------------------------------------------------
!     Apply boundary deformation restrictions
! ------------------------------------------------------------------------------  
  
      patchcheckcounter = 0
        
! --- 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%movePatchDir
      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%movePatchDir
      END DO ! iPatch                   
    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
! Reduction of rhs 
! END TO DO 
                        
! ==============================================================================
!   Update displacements
! ==============================================================================

    DO ireg = 1,global%nRegionsLocal
      pgrid => regions(ireg)%grid
      pdisp => pgrid%disp
      prhs  => pgrid%rhs    
    
      DO iv = 1,pgrid%nVertTot
        term = 1.0_rfreal/REAL(pGrid%gmVertWght(iv),RFREAL)     
      
        pdisp(xcoord,iv) = pdisp(xcoord,iv) + term*prhs(xcoord,iv)
        pdisp(ycoord,iv) = pdisp(ycoord,iv) + term*prhs(ycoord,iv)
        pdisp(zcoord,iv) = pdisp(zcoord,iv) + term*prhs(zcoord,iv)
      END DO ! iv      
    END DO ! iReg

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

! ******************************************************************************
! Update coordinates
! ******************************************************************************

  IF ( niter > 0 ) THEN 
    DO ireg = 1,global%nRegionsLocal 
      pgrid => regions(ireg)%grid
      pdisp => pgrid%disp
      pxyz  => pgrid%xyz
    
      DO iv = 1,pgrid%nVertTot
        pxyz(xcoord,iv) = pxyz(xcoord,iv) + pdisp(xcoord,iv)
        pxyz(ycoord,iv) = pxyz(ycoord,iv) + pdisp(ycoord,iv)
        pxyz(zcoord,iv) = pxyz(zcoord,iv) + pdisp(zcoord,iv)            
      END DO ! iv
    END DO ! iReg
  END IF ! nIter

! ******************************************************************************
! Deallocate temporary arrays
! ******************************************************************************  
  
  IF ( niter > 0 ) THEN 
    DO ireg = 1,global%nRegionsLocal 
      pgrid => regions(ireg)%grid
          
      DEALLOCATE(pgrid%gmEdgeWght,stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_deallocate,__line__,'pGrid%gmEdgeWght')
      END IF ! global%error        
          
      DEALLOCATE(pgrid%gmVertWght,stat=errorflag)
      global%error = errorflag
      IF ( global%error /= err_none ) THEN 
        CALL errorstop(global,err_deallocate,__line__,'pGrid%gmVertWght')
      END IF ! global%error               
    END DO ! iReg  
  END IF ! nIter

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

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

  CALL deregisterfunction(global)

END SUBROUTINE rflu_movegriddisp

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_MoveGridDisp.F90,v $
! Revision 1.11  2008/12/06 08:44:30  mtcampbe
! Updated license.
!
! Revision 1.10  2008/11/19 22:17:42  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.9  2005/06/09 20:22:58  haselbac
! Replaced movePatch by movePatchDir
!
! Revision 1.8  2005/04/15 15:07:20  haselbac
! Converted to MPI, grid motion can only be used for serial runs
!
! Revision 1.7  2004/09/17 19:59:58  haselbac
! Added IF statements to prevent operations for zero iterations
!
! Revision 1.6  2004/04/14 03:57:11  haselbac
! Bug fix: Now get correct boundary deformation for zero smoothing iterations
!
! Revision 1.5  2003/08/20 02:09:58  haselbac
! Changed verbosity conditions to reduce solver output in GENx runs
!
! Revision 1.4  2003/07/09 22:36:53  haselbac
! Added IF statement to avoid probs with single-proc charm jobs
!
! Revision 1.3  2003/05/01 14:12:10  haselbac
! Added logic to deal with non-moving and non-smoothed patches
!
! Revision 1.2  2003/04/17 00:16:19  haselbac
! Bug fix: Weight should be inverse sqrt, little effect
!
! Revision 1.1  2003/03/31 16:05:39  haselbac
! Initial revision
!
! ******************************************************************************