PLAG_GetCellIndices.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: Search algorithm for cell indices.
!
! Description: none.
!
! Input: region = current region.
!
! Output: region%levels%plag%aiv
!         region%levels%plag%aivOld
!
! Notes: none.
!
!******************************************************************************
!
! $Id: PLAG_GetCellIndices.F90,v 1.4 2009/10/26 00:19:32 mtcampbe Exp $
!
! Copyright: (c) 2002 by the University of Illinois
!
!******************************************************************************

SUBROUTINE plag_getcellindices( region, iReg )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modpartlag, ONLY    : t_plag, t_plag_input
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
                            rflo_getnodeoffset
  USE plag_modinterfaces, ONLY : plag_incelltest, plag_incelltestrobust
  USE moderror  
  USE modparameters
  USE plag_modparameters
  IMPLICIT NONE

#include "Indexing.h"

! ... parameters
  TYPE(t_region) :: region

  INTEGER :: ireg

! ... loop variables
  INTEGER :: i, icelllev, ilevels, ipcls, j, k

! ... local variables
  CHARACTER(CHRLEN) :: rcsidentstring

  INTEGER :: icoff, ijcoff, ijkc, inoff, ijnoff, ilev,                        &
             ijknr, ijknri, ijknrj,ijknrk, ihigh, ilow, iskip, ipcbeg,ipcend, &
             jpcbeg,jpcend, jhigh, jlow, jskip,                               &
             kpcbeg,kpcend, khigh, klow, kskip,                               &   
             ncont, ncelllev, ncelllevmax, nlevels, npcls 
  INTEGER,          DIMENSION(4)   :: indexcurr, indexnew, indexsearch
  INTEGER, POINTER, DIMENSION(:,:) :: paiv, paivold
  
  LOGICAL :: celllocate, celllocaterobust
  
  REAL(RFREAL) :: massl,taulr,diaml

  REAL(RFREAL),          DIMENSION(3)   :: posplag
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pcv, pcvold
  REAL(RFREAL), POINTER, DIMENSION(:,:) :: pdv
  
  TYPE(t_plag),   POINTER :: pplag  
  TYPE(t_global), POINTER :: global
  TYPE(t_region), POINTER :: pregion  
  
!******************************************************************************

  rcsidentstring = '$RCSfile: PLAG_GetCellIndices.F90,v $ $Revision: 1.4 $'

  global => region%global
  
  CALL registerfunction( global, 'PLAG_getCellIndices',&
  'PLAG_GetCellIndices.F90' )

! Get dimensions --------------------------------------------------------------

  ilev  = region%currLevel

  npcls = region%levels(ilev)%plag%nPcls 
  ncelllevmax = 1
  nlevels     = 8
  
  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
  
! Set pointers ----------------------------------------------------------------
   
  pplag   => region%levels(ilev)%plag

  paiv    => pplag%aiv 
  paivold => pplag%aivOld
  pcv     => pplag%cv
  pcvold  => pplag%cvOld
  pdv     => pplag%dv
 
! Loop over Lagrangian particles ----------------------------------------------

  DO ipcls = 1, npcls      

! - Load current index --------------------------------------------------------
    
    indexcurr(1) = paivold(aiv_plag_icells,ipcls) 
    indexcurr(2) = paivold(aiv_plag_indexi,ipcls) 
    indexcurr(3) = paivold(aiv_plag_indexj,ipcls) 
    indexcurr(4) = paivold(aiv_plag_indexk,ipcls) 

! - Load particle positions ---------------------------------------------------

    posplag(xcoord) = pcv(cv_plag_xpos,ipcls)
    posplag(ycoord) = pcv(cv_plag_ypos,ipcls)
    posplag(zcoord) = pcv(cv_plag_zpos,ipcls)
 
! - Set particle status -------------------------------------------------------

       paiv(aiv_plag_status,ipcls) = plag_status_keep
    paivold(aiv_plag_status,ipcls) = plag_status_keep

! - Search current and surrounding cells --------------------------------------
    
    ncelllev = 1
    
    DO WHILE(ncelllev <= ncelllevmax)
    
      DO ilevels = 1, nlevels 
      
        SELECT CASE(ilevels)
        
! -- current cell -------------------------------------------------------------
          CASE(1)
            ilow  = indexcurr(2)
            ihigh = indexcurr(2)
            iskip = 1
            
            jlow  = indexcurr(3)
            jhigh = indexcurr(3)
            jskip = 1 
            
            klow  = indexcurr(4)
            khigh = indexcurr(4)
            kskip = 1

! -- i-cell shift -------------------------------------------------------------
          CASE(2)
            ilow  = indexcurr(2)-ncelllev
            ihigh = indexcurr(2)+ncelllev
            iskip = ncelllev+1
            
            jlow  = indexcurr(3)
            jhigh = indexcurr(3)
            jskip = 1 
            
            klow  = indexcurr(4)
            khigh = indexcurr(4)
            kskip = 1

! -- j-cell shift -------------------------------------------------------------
          CASE(3)
            ilow  = indexcurr(2)
            ihigh = indexcurr(2)
            iskip = 1
            
            jlow  = indexcurr(3)-ncelllev
            jhigh = indexcurr(3)+ncelllev
            jskip = ncelllev+1
            
            klow  = indexcurr(4)
            khigh = indexcurr(4)
            kskip = 1
 
! -- k-cell shift -------------------------------------------------------------
          CASE(4)
            ilow  = indexcurr(2)
            ihigh = indexcurr(2)
            iskip = 1
            
            jlow  = indexcurr(3)
            jhigh = indexcurr(3)
            jskip = 1
            
            klow  = indexcurr(4)-ncelllev
            khigh = indexcurr(4)+ncelllev
            kskip = ncelllev+1
                           
! -- ij edge cell shift -------------------------------------------------------------
          CASE(5)
            ilow  = indexcurr(2)-ncelllev
            ihigh = indexcurr(2)+ncelllev
            iskip = ncelllev+1
    
            jlow  = indexcurr(3)-ncelllev
            jhigh = indexcurr(3)+ncelllev
            jskip = ncelllev+1
     
            klow  = indexcurr(4)
            khigh = indexcurr(4)
            kskip = 1
    
! -- ik edge cell shift -------------------------------------------------------
          CASE(6)
            ilow  = indexcurr(2)-ncelllev
            ihigh = indexcurr(2)+ncelllev
            iskip = ncelllev+1
            
            jlow  = indexcurr(3)
            jhigh = indexcurr(3)
            jskip = 1
            
            klow  = indexcurr(4)-ncelllev
            khigh = indexcurr(4)+ncelllev
            kskip = ncelllev+1

! -- jk edge cell shift -------------------------------------------------------
          CASE(7) 
            ilow  = indexcurr(2)
            ihigh = indexcurr(2)
            iskip = 1
                           
            jlow  = indexcurr(3)-ncelllev
            jhigh = indexcurr(3)+ncelllev
            jskip = ncelllev+1
     
            klow  = indexcurr(4)-ncelllev
            khigh = indexcurr(4)+ncelllev
            kskip = ncelllev+1

! -- ijk corner cell shift ----------------------------------------------------  
          CASE(8) 
            ilow  = indexcurr(2)-ncelllev
            ihigh = indexcurr(2)+ncelllev
            iskip = ncelllev+1
     
            jlow  = indexcurr(3)-ncelllev
            jhigh = indexcurr(3)+ncelllev
            jskip = ncelllev+1
    
            klow  = indexcurr(4)-ncelllev
            khigh = indexcurr(4)+ncelllev
            kskip = ncelllev+1
            
        END SELECT !iLevels 

! -- Perform search -----------------------------------------------------------
           
        DO k=klow,khigh,kskip
        DO j=jlow,jhigh,jskip
        DO i=ilow,ihigh,iskip
          ijkc   = indijk(i,  j  ,k  ,icoff,ijcoff)
          ijknr  = indijk(i,  j  ,k  ,inoff,ijnoff)              
          ijknri = indijk(i+1,j  ,k  ,inoff,ijnoff)
          ijknrj = indijk(i,j+1  ,k  ,inoff,ijnoff)
          ijknrk = indijk(i  ,j  ,k+1,inoff,ijnoff)

          indexsearch(1) = ijkc
          indexsearch(2) = i
          indexsearch(3) = j
          indexsearch(4) = k
  
          CALL plag_incelltest(region, posplag, indexsearch, &
                               ijknr,ijknri,ijknrj,ijknrk,   &
                               indexnew,celllocate)

          IF ( celllocate .EQV. .true. ) goto 999
        
        END DO ! i 
        END DO ! j
        END DO ! k

      END DO ! iLevels  
        
      ncelllev = ncelllev+1

    END DO ! nCellLev

! - Trap error if unable to locate --------------------------------------------
     
    IF ( celllocate .EQV. .false. ) THEN       
      massl = sum( pplag%cv(pplag%cvPlagMass(:),ipcls) )
      diaml = pdv(dv_plag_diam,ipcls)

      IF ( diaml > 1.0e-14_rfreal ) THEN
        taulr = 3.0_rfreal*global%pi*pplag%tv(tv_plag_muelmixt,ipcls)*diaml/massl
      ELSE
        taulr = 0.0_rfreal
      ENDIF

      WRITE(stdout,'(A)') &
       '##### Rocpart Warning: PLAG_inCellTest Unable to Locate Cell #####'   
      WRITE(stdout,1010) 'Default Search Failed to Locate Cell for Particle',&
       global%currentTime, ireg, ipcls,&
       paivold(aiv_plag_pidini,ipcls),          &
       paivold(aiv_plag_regini,ipcls),          &
       paivold(aiv_plag_icells,ipcls),          &
       paivold(aiv_plag_indexi,ipcls),          &
       paivold(aiv_plag_indexj,ipcls),          &
       paivold(aiv_plag_indexk,ipcls),          &
       pcvold(cv_plag_xpos:cv_plag_zpos,ipcls), &
       posplag(1:3),                            &
       pdv(dv_plag_uvel:dv_plag_wvel,ipcls),    &
       pdv(dv_plag_uvelmixt:dv_plag_wvelmixt,ipcls),&
       pdv(dv_plag_diam,ipcls),massl,taulr
    
! -- Invoke robust cell search for particles potentially located --------------
!      in dummy cells ---------------------------------------------------------

      CALL plag_incelltestrobust( region, posplag,indexcurr, &
                                  indexnew,celllocaterobust  )

! TEMPORARY
! Delete Particle from data structure failing robust cell search
!  due to clobbered geometry datastructure. Relevant for complex geometries.
!
      WRITE(stdout,'(A)') &
       '##### Rocpart Warning: Deleting Particle Following Robust Cell Search #####'
         paiv(aiv_plag_status,ipcls) = plag_status_delete
      paivold(aiv_plag_status,ipcls) = plag_status_delete
! END TEMPORARY

      IF ( celllocaterobust .EQV. .true. ) goto 999
      
        WRITE(stdout,'(A)') &
         '##### Rocpart Error: PLAG_inCellTestRobust Unable to Locate Cell with Robust Search #####'   
        WRITE(stdout,1010) 'Unable to Locate Cell for Particle',&
        global%currentTime, ireg, ipcls,&
        paivold(aiv_plag_pidini,ipcls),          &
        paivold(aiv_plag_regini,ipcls),          &
        paivold(aiv_plag_icells,ipcls),          &
        paivold(aiv_plag_indexi,ipcls),          &
        paivold(aiv_plag_indexj,ipcls),          &
        paivold(aiv_plag_indexk,ipcls),          &
        pcvold(cv_plag_xpos:cv_plag_zpos,ipcls), &
        posplag(1:3),                            &
        pdv(dv_plag_uvel:dv_plag_wvel,ipcls),    &
        pdv(dv_plag_uvelmixt:dv_plag_wvelmixt,ipcls)
      
!        CALL ErrorStop( global,ERR_PLAG_CELLINDEX,__LINE__ )

! TEMPORARY
! Delete Particle from data structure for testing
         paiv(aiv_plag_status,ipcls) = plag_status_delete 
      paivold(aiv_plag_status,ipcls) = plag_status_delete 
! END TEMPORARY

    END IF ! cellLocate 

! - Load new cell indices -----------------------------------------------------
         
999 CONTINUE

!  Dont load aiv into aivOld till end of RK-stage since 
!   search algorithm is modulated by RK time-stepping that
!   is based on cvOld
!
!    pAivOld(AIV_PLAG_ICELLS,iPcls) = pAiv(AIV_PLAG_ICELLS,iPcls)
!    pAivOld(AIV_PLAG_INDEXI,iPcls) = pAiv(AIV_PLAG_INDEXI,iPcls)
!    pAivOld(AIV_PLAG_INDEXJ,iPcls) = pAiv(AIV_PLAG_INDEXJ,iPcls)
!    pAivOld(AIV_PLAG_INDEXK,iPcls) = pAiv(AIV_PLAG_INDEXK,iPcls)

    paiv(aiv_plag_icells,ipcls) = indexnew(1)
    paiv(aiv_plag_indexi,ipcls) = indexnew(2) 
    paiv(aiv_plag_indexj,ipcls) = indexnew(3) 
    paiv(aiv_plag_indexk,ipcls) = indexnew(4)
  
  END DO  ! iPcls 
   
! finalize --------------------------------------------------------------------

  CALL deregisterfunction( global )

1010 FORMAT(a,'Time = ', 1pe12.5, ', iReg = ', i5, ', iPcls = ', i8, &
            ', aivOld = ',6i5,', posOld = ',3(1pe15.7),&
            ', posCurr = ', 3(1pe15.7),', pVel = ', 3(1pe15.7), &
            ', mVel = ', 3(1pe15.7),&
            ', diam = ', 1pe15.7,', mass = ',1pe15.7,' , tauL = ', 1pe15.7)

END SUBROUTINE plag_getcellindices
!******************************************************************************
!
! RCS Revision history:
!
! $Log: PLAG_GetCellIndices.F90,v $
! Revision 1.4  2009/10/26 00:19:32  mtcampbe
! Updates for completion of NATIVE_MP_IO
!
! Revision 1.3  2008/12/06 08:44:33  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:46  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/12/01 20:57:33  fnajjar
! Initial revision after changing case
!
! Revision 1.9  2004/08/05 20:10:41  fnajjar
! Included deletion flag for aiv after robust cell search algorithm
!
! Revision 1.8  2004/07/01 14:15:54  fnajjar
! Added variables to IO when particle search fails
!
! Revision 1.7  2004/07/01 14:13:38  fnajjar
! Commented out aivOld loading from aiv since it should be done at initial RK-stage
!
! Revision 1.6  2004/04/09 23:10:20  fnajjar
! Added robust kernel for incell testing
!
! Revision 1.5  2003/05/15 02:57:05  jblazek
! Inlined index function.
!
! Revision 1.4  2003/04/18 22:49:04  fnajjar
! Bug fix to move IF statement of cellLocation inside ijk DO-loop
!
! Revision 1.3  2003/01/16 20:43:27  f-najjar
! Include iReg in calling sequence
!
! Revision 1.2  2003/01/16 20:15:11  f-najjar
! Removed iRegionGlobal
!
! Revision 1.1  2002/10/25 14:15:43  f-najjar
! Initial Import of Rocpart
!
!
!******************************************************************************