RFLO_CheckValidity.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: Check validity of variables
!
! Description: It detects NaN solution and/or negative pressure or density.
!
! Input:  region = Region data
!
! Output: none.
!
! Notes: none.
!
!******************************************************************************
!
! $Id: RFLO_CheckValidity.F90,v 1.6 2008/12/06 08:44:06 mtcampbe Exp $
!
! Copyright: (c) 2003 by the University of Illinois
!
!******************************************************************************

SUBROUTINE rflo_checkvalidity(region)

  USE moddatatypes
  USE moderror
  USE modglobal, ONLY: t_global
  USE moddatastruct, ONLY: t_region
  USE modbndpatch, ONLY: t_patch
  USE modparameters
  USE modmpi
  USE modtools, ONLY: isnan
  USE modinterfaces, ONLY: mixtperf_g_cpr, mixtperf_p_deogvm2, mixtperf_r_m, &
                           mixtperf_t_dpr, &
                           rflo_getdimensphys, rflo_getcelloffset
  
  IMPLICIT NONE
#include "Indexing.h"

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

! ... loop variables
  INTEGER :: i, j, k, ipatch

! ... local variables
  TYPE(t_global), POINTER :: global
  TYPE(t_patch), POINTER  :: patches(:)

  INTEGER, PARAMETER :: max_invalid_locs = 10
  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ilev, icoff, ijcoff, ijkc
  INTEGER :: indcp,indmol,nlocs
#ifdef MPI
  INTEGER :: mpierrcode
#endif
  REAL(RFREAL) :: eo,gamma,p,rgas,rho,rrho,t,u,v,vm2,w,rmin,pmin
  REAL(RFREAL), DIMENSION(:,:), POINTER :: cv,gv
  LOGICAL :: foundnan

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

  global => region%global

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

  nlocs = 0

! get dimensions and pointers -------------------------------------------------

  ilev = region%currLevel

  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )

  cv      => region%levels(ilev)%mixt%cv
  gv      => region%levels(ilev)%mixt%gv
  patches => region%levels(ilev)%patches

  indcp  = region%levels(ilev)%mixt%indCp
  indmol = region%levels(ilev)%mixt%indMol 

! loop over cells and check for positivity ------------------------------------

  rmin = 1.e+13_rfreal
  pmin = 1.e+13_rfreal
  foundnan = .false.

  DO k=kpcbeg,kpcend
    DO j=jpcbeg,jpcend
      DO i=ipcbeg,ipcend
        ijkc = indijk(i,j,k,icoff,ijcoff)

        rho  = cv(cv_mixt_dens,ijkc)
        rrho = 1.0_rfreal/rho
        u    = rrho*cv(cv_mixt_xmom,ijkc)
        v    = rrho*cv(cv_mixt_ymom,ijkc)
        w    = rrho*cv(cv_mixt_zmom,ijkc)        
        eo   = rrho*cv(cv_mixt_ener,ijkc)

        rgas  = mixtperf_r_m(gv(gv_mixt_mol,ijkc*indmol))
        gamma = mixtperf_g_cpr(gv(gv_mixt_cp,ijkc*indcp),rgas)
        vm2   = u*u + v*v + w*w

        p = mixtperf_p_deogvm2(rho,eo,gamma,vm2)

        rmin  = min( rmin,rho )
        pmin  = min( pmin,p )

        foundnan = isnan(rho)
        IF (foundnan) goto 888

      ENDDO      ! i
    ENDDO        ! j
  ENDDO          ! k

888 CONTINUE

  IF (foundnan .OR. rmin < 0._rfreal .OR. pmin < 0._rfreal) THEN

    DO k=kpcbeg,kpcend
     DO j=jpcbeg,jpcend
      DO i=ipcbeg,ipcend
        ijkc = indijk(i,j,k,icoff,ijcoff)

        rho  = cv(cv_mixt_dens,ijkc)
        rrho = 1.0_rfreal/rho
        u    = rrho*cv(cv_mixt_xmom,ijkc)
        v    = rrho*cv(cv_mixt_ymom,ijkc)
        w    = rrho*cv(cv_mixt_zmom,ijkc)        
        eo   = rrho*cv(cv_mixt_ener,ijkc)

        rgas  = mixtperf_r_m(gv(gv_mixt_mol,ijkc*indmol))
        gamma = mixtperf_g_cpr(gv(gv_mixt_cp,ijkc*indcp),rgas)
        vm2   = u*u + v*v + w*w

        p = mixtperf_p_deogvm2(rho,eo,gamma,vm2)
        t = mixtperf_t_dpr(rho,p,rgas)

        IF ( (isnan(rho) .EQV. .true.) .OR. & 
             (isnan(u)   .EQV. .true.) .OR. &
             (isnan(v)   .EQV. .true.) .OR. & 
             (isnan(w)   .EQV. .true.) .OR. & 
             (isnan(p)   .EQV. .true.) .OR. &
             (isnan(t)   .EQV. .true.) .OR. &
             (      rho  <  0._rfreal) .OR. &
             (      p    <  0._rfreal) ) THEN
          nlocs = nlocs + 1   

          IF ( nlocs == 1 ) THEN 
            WRITE(stdout,'(A,1X,A,1X,I9)') solver_name, & 
                  'Invalid variables detected!'
              
            IF ( global%flowType == flow_unsteady ) THEN 
              WRITE(stdout,'(A,3X,A,1X,1PE12.5)') solver_name, &
                                                  'Current time:', &
                                                  global%currentTime              
            ELSE 
              WRITE(stdout,'(A,3X,A,1X,I6)') solver_name, &
                                             'Current iteration number:', &
                                             global%currentIter           
            END IF ! global%flowType                 
                                            
            WRITE(stdout,1000) solver_name,'Region:',region%iRegionGlobal, & 
                        ', bc-types:', (patches(ipatch)%bcType, &
                          ipatch = 1,region%nPatches)
                                  
            WRITE(stdout,'(A,6X,A,7(1X,A))') solver_name,'#', &
                                             '   Density   ', &
                                             '  x-velocity ', &
                                             '  y-velocity ', &
                                             '  z-velocity ', &
                                             '   Pressure  ', &
                                             ' Temperature ', &
                                             '   i   ,   j   ,   k'       
          END IF ! nLocs

          IF ( nlocs <= max_invalid_locs ) THEN 
            WRITE(stdout,'(A,4X,I3,6(1X,E13.6),2X,3I7)') solver_name,nlocs, & 
                                                  rho,u,v,w,p,t,i,j,k
          END IF ! nLocs
        END IF   ! dv    
      ENDDO      ! i
     ENDDO       ! j
    ENDDO        ! k

#ifdef MPI
    CALL mpi_abort( global%mpiComm,mpierrcode,global%mpierr )
#endif
    stop
  ENDIF          ! foundNan...

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

1000 FORMAT( a,3x,a,i5,a,60i5 )

  CALL deregisterfunction( global )

END SUBROUTINE rflo_checkvalidity

!******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLO_CheckValidity.F90,v $
! Revision 1.6  2008/12/06 08:44:06  mtcampbe
! Updated license.
!
! Revision 1.5  2008/11/19 22:17:20  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.4  2005/05/12 21:31:14  wasistho
! write i,j,k in the same line as previous vars
!
! Revision 1.3  2005/05/12 20:53:51  wasistho
! added bctypes in error msg
!
! Revision 1.2  2005/01/11 00:26:15  wasistho
! changed mpi_finalize to mpi_abort
!
! Revision 1.1  2004/11/29 21:25:16  wasistho
! lower to upper case
!
! Revision 1.6  2004/08/24 01:02:23  wasistho
! search nan per cell i.o. the sum value
!
! Revision 1.5  2004/08/16 17:05:52  wasistho
! moved MPI finalize within IF statement
!
! Revision 1.4  2004/08/04 00:29:35  wasistho
! speedup check validity
!
! Revision 1.3  2004/07/26 20:03:13  wasistho
! added i,j,k locations
!
! Revision 1.2  2004/07/26 19:30:44  wasistho
! changed POINTER to TARGET for region parameter
!
! Revision 1.1  2004/07/26 19:10:06  wasistho
! initial import RFLO_CheckValidity
!
!
!******************************************************************************