PLAG_WallBounce.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: perform elastic bounce of particles
!          for slip wall, non-slip wall, injection 
!          and symmetry boundary conditions. 
!
! Description: none.
!
! Input: region = current region.
!
! Output: region%levels%plag%cv = plag data.
!
! Notes:
!
!   The cell index for the particle has not yet been updated to its position
!
!******************************************************************************
!
! $Id: PLAG_WallBounce.F90,v 1.3 2008/12/06 08:44:36 mtcampbe Exp $
!
! Copyright: (c) 2002 by the University of Illinois
!
!******************************************************************************

SUBROUTINE plag_wallbounce( region )

  USE moddatatypes
  USE modpartlag, ONLY    : t_plag, t_plag_input
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global  
  USE modinterfaces, ONLY : rflo_getpatchindices, rflo_getpatchdirection, &
                            rflo_getnodeoffset,   rflo_getdimensphys
  USE moderror
  USE modparameters
  USE plag_modparameters
  IMPLICIT NONE

#include "Indexing.h"

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

! ... loop variables
  INTEGER :: ipatch, ipcls

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

  INTEGER :: bctype, ilev, lbound, npatches, npcls
  INTEGER :: inoff, ijnoff, ijknpatch
  INTEGER :: icplag, ipcbeg, ipcend, ibeg, iend, idir
  INTEGER :: jcplag, jpcbeg, jpcend, jbeg, jend, jdir
  INTEGER :: kcplag, kpcbeg, kpcend, kbeg, kend, kdir

  INTEGER, POINTER, DIMENSION(:,:) :: paiv

  LOGICAL :: lboundskip(6)

  REAL(RFREAL)               :: dpgrid, dpgridold, dpmome, dpmomeold, dfac, sgn
  REAL(RFREAL)               :: dpmomerhssum,dpposrhssum
  REAL(RFREAL), DIMENSION(3) :: diffpos, facecentroid, momeplag, posplag, sface
  REAL(RFREAL), DIMENSION(3) :: diffposold, momeplagold, posplagold
  REAL(RFREAL), DIMENSION(3) :: momerhssumplag, posrhssumplag
  REAL(RFREAL), POINTER, DIMENSION(:,:)   :: pcv, pcvold, prhssum, psi, psj, psk
  REAL(RFREAL), POINTER, DIMENSION(:,:,:) :: pfc

  TYPE(t_patch),  POINTER :: ppatch
  TYPE(t_plag),   POINTER :: pplag 
  TYPE(t_global), POINTER :: global

!******************************************************************************

  rcsidentstring = '$RCSfile: PLAG_WallBounce.F90,v $ $Revision: 1.3 $'

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

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

  ilev  = region%currLevel
  npcls = region%levels(ilev)%plag%nPcls 
  
  IF (npcls == 0) goto 999 ! exit if no particles

  npatches = region%nPatches

  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )

! Set pointers ----------------------------------------------------------------

  pplag  => region%levels(ilev)%plag 
  paiv   => pplag%aiv  
  pcv    => pplag%cv
  pcvold => pplag%cvOld
  prhssum => pplag%rhsSum

  psi    => pplag%si
  psj    => pplag%sj
  psk    => pplag%sk

  pfc    => pplag%fc 

! Loop over particles ---------------------------------------------------------

  DO ipcls = 1, npcls

    icplag = paiv(aiv_plag_indexi,ipcls)
    jcplag = paiv(aiv_plag_indexj,ipcls)
    kcplag = paiv(aiv_plag_indexk,ipcls)

! - Cycle if particle is not adjacent to a boundary ---------------------------

    IF (ipcbeg < icplag .AND. icplag < ipcend .AND. &
        jpcbeg < jcplag .AND. jcplag < jpcend .AND. &
        kpcbeg < kcplag .AND. kcplag < kpcend ) cycle

    posplag(1:3) = pcv(cv_plag_xpos:cv_plag_zpos,ipcls)

! - Set lboundSkip(:) to .TRUE. for values of lbound adjacent to a boundary ---

    lboundskip(1:6) = .true.

    IF (icplag <= ipcbeg) lboundskip(1) = .false. ! to include ipcbeg = ipcend
    IF (icplag >= ipcend) lboundskip(2) = .false. ! case, ELSE IF not used here

    IF (jcplag <= jpcbeg) lboundskip(3) = .false.
    IF (jcplag >= jpcend) lboundskip(4) = .false.

    IF (kcplag <= kpcbeg) lboundskip(5) = .false.
    IF (kcplag >= kpcend) lboundskip(6) = .false.

! - Loop over patches ---------------------------------------------------------

    DO ipatch=1,npatches

      ppatch => region%levels(ilev)%patches(ipatch)
      lbound =  ppatch%lbound
      IF (lboundskip(lbound)) cycle

      bctype = ppatch%bcType

! --- Select specific boundary condition types --------------------------------

      IF ((bctype>=bc_slipwall   .AND. bctype<=bc_slipwall  +bc_range) .OR. &
          (bctype>=bc_noslipwall .AND. bctype<=bc_noslipwall+bc_range) .OR. &
          (bctype>=bc_injection  .AND. bctype<=bc_injection +bc_range) .OR. &
          (bctype>=bc_symmetry   .AND. bctype<=bc_symmetry  +bc_range)) THEN 

! ----- Check if particle cell is within (physical cells of) patch ------------

        CALL rflo_getpatchindices( region,ppatch,ilev, &
                                   ibeg,iend,jbeg,jend,kbeg,kend )

        IF ( ibeg <= icplag .AND. icplag <= iend .AND. &
             jbeg <= jcplag .AND. jcplag <= jend .AND. &
             kbeg <= kcplag .AND. kcplag <= kend ) THEN

! ------- Select correct face vector and make it point inwards ----------------

          IF (lbound==1 .OR. lbound==3 .OR. lbound==5) THEN
            sgn   = -1.0_rfreal
            ijknpatch = indijk(icplag,jcplag,kcplag,inoff,ijnoff)
          ELSE
            sgn   = +1.0_rfreal
            CALL rflo_getpatchdirection( ppatch,idir,jdir,kdir )
            ijknpatch = indijk(icplag-idir,jcplag-jdir,kcplag-kdir,inoff,ijnoff)
          ENDIF ! lbound

          SELECT CASE (lbound)

            CASE(1,2)
              sface(1:3)        = sgn*psi(xcoord:zcoord,       ijknpatch)
              facecentroid(1:3) =     pfc(xcoord:zcoord,icoord,ijknpatch)

            CASE(3,4)
              sface(1:3)        = sgn*psj(xcoord:zcoord,       ijknpatch)
              facecentroid(1:3) =     pfc(xcoord:zcoord,jcoord,ijknpatch)

            CASE(5,6)
              sface(1:3)        = sgn*psk(xcoord:zcoord,       ijknpatch)
              facecentroid(1:3) =     pfc(xcoord:zcoord,kcoord,ijknpatch)

          END SELECT ! lbound

          diffpos(1:3) = posplag(1:3)-facecentroid(1:3)
          dpgrid       = dot_product(sface,diffpos)

! ------- Apply reflection to particle exiting computational domain -----------

          IF ( dpgrid < 0.0_rfreal ) THEN

            momeplag(1:3) = pcv(cv_plag_xmom:cv_plag_zmom,ipcls)
            dpmome = dot_product(sface,momeplag)
            dfac   = -2.0_rfreal/dot_product(sface,sface)
            pcv(cv_plag_xpos:cv_plag_zpos,ipcls) = posplag(1:3)  + &
                                                   dpgrid*dfac*sface(1:3)
            pcv(cv_plag_xmom:cv_plag_zmom,ipcls) = momeplag(1:3) + &
                                                   dpmome*dfac*sface(1:3)

! --------- Also need to reflect cvOld values:  the particle must seem like ---
! --------- it is passing straight through the boundary from the other side ---

            posplagold(1:3)  = pcvold(cv_plag_xpos:cv_plag_zpos,ipcls)
            diffposold(1:3)  = posplagold(1:3)-facecentroid(1:3)
            dpgridold        = dot_product(sface,diffposold)

            momeplagold(1:3) = pcvold(cv_plag_xmom:cv_plag_zmom,ipcls)
            dpmomeold        = dot_product(sface,momeplagold)

            pcvold(cv_plag_xpos:cv_plag_zpos,ipcls) = posplagold(1:3)  + &
                                                      dpgridold*dfac*sface(1:3)
            pcvold(cv_plag_xmom:cv_plag_zmom,ipcls) = momeplagold(1:3) + &
                                                      dpmomeold*dfac*sface(1:3)

! --------- Apply reflectivity to rhsSum of positions and momenta for consistency 

            posrhssumplag(1:3)  = prhssum(cv_plag_xpos:cv_plag_zpos,ipcls)
            dpposrhssum         = dot_product(sface,posrhssumplag)

            momerhssumplag(1:3) = prhssum(cv_plag_xmom:cv_plag_zmom,ipcls)
            dpmomerhssum        = dot_product(sface,momerhssumplag)

            prhssum(cv_plag_xpos:cv_plag_zpos,ipcls) = posrhssumplag(1:3)  &
                                                     + dpposrhssum*dfac*sface(1:3)
            prhssum(cv_plag_xmom:cv_plag_zmom,ipcls) = momerhssumplag(1:3) &
                                                     + dpmomerhssum*dfac*sface(1:3)

! --------- Note: no effect on energy because bounce is purely elastic --------

          ENDIF ! dpGrid

        ENDIF   ! iCPlag, jCPlag, kCPlag
      ENDIF     ! bcType
    ENDDO       ! iPatch
  ENDDO         ! iPcls

! finalize --------------------------------------------------------------------

999  CONTINUE
  CALL deregisterfunction( global )

END SUBROUTINE plag_wallbounce

!******************************************************************************
!
! RCS Revision history:
!
! $Log: PLAG_WallBounce.F90,v $
! Revision 1.3  2008/12/06 08:44:36  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:48  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/12/01 20:58:21  fnajjar
! Initial revision after changing case
!
! Revision 1.8  2004/06/17 14:32:16  fnajjar
! Applied reflectivity to rhsSum vectors of positions and momenta
!
! Revision 1.7  2004/03/01 16:38:00  jferry
! added reflection of cvOld quantities
!
! Revision 1.6  2003/11/03 21:21:51  fnajjar
! Changed definition of face vectors pointing to PLAG datastructure
!
! Revision 1.5  2003/05/15 02:57:05  jblazek
! Inlined index function.
!
! Revision 1.4  2003/05/01 22:58:30  jferry
! overhauled structure in order to optimize performance
!
! Revision 1.3  2003/04/16 23:23:25  fnajjar
! Bug fix for appropriate normals of face vectors
!
! Revision 1.2  2003/01/16 20:27:08  f-najjar
! Removed iRegionGlobal
!
! Revision 1.1  2002/10/25 14:20:32  f-najjar
! Initial Import of Rocpart
!
!******************************************************************************