PEUL_EnforcePositivity.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 if there are negative values of concentration, and, when
!          indicated, report and possibly fix these values.
!
! Description: none.
!
! Input: region = data of current region.
!
! Output: region%levels%peul%cv
!
! Notes: none.
!
!******************************************************************************
!
! $Id: PEUL_EnforcePositivity.F90,v 1.3 2008/12/06 08:44:39 mtcampbe Exp $
!
! Copyright: (c) 2004 by the University of Illinois
!
!******************************************************************************

SUBROUTINE peul_enforcepositivity( region )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal,     ONLY : t_global
  USE moderror
  USE modparameters
  USE peul_modparameters

  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset
  IMPLICIT NONE

#include "Indexing.h"

! ... parameters
  TYPE(t_region), INTENT(INOUT) :: region

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

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

  INTEGER :: ipcbeg,ipcend,jpcbeg,jpcend,kpcbeg,kpcend,nptypes
  INTEGER :: ilev,icoff,ijcoff,ijkc0
  INTEGER :: ncnt1,ncnt2,iloc1,jloc1,kloc1,iloc2,jloc2,kloc2

  LOGICAL :: clip,firsttype

  REAL(RFREAL) :: minc,maxc,maxconc
  REAL(RFREAL), POINTER :: cv(:,:)

  TYPE(t_global), POINTER :: global

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

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

  global => region%global

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

! begin -----------------------------------------------------------------------

! 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)%peul%cv

  nptypes = region%peulInput%nPtypes

  firsttype = .true.

  DO ipt=1,nptypes

    maxconc = region%peulInput%ptypes(ipt)%maxConc

    minc =  1.e20_rfreal
    maxc = -1.e20_rfreal
    ncnt1 = 0
    ncnt2 = 0
    iloc1 = crazy_value_int
    jloc1 = crazy_value_int
    kloc1 = crazy_value_int
    iloc2 = crazy_value_int
    jloc2 = crazy_value_int
    kloc2 = crazy_value_int

    clip = .false.

    SELECT CASE (region%peulInput%ptypes(ipt)%clipModel)

    CASE (peul_clip_model_used)
      clip = .true.

    END SELECT ! clipModel

    IF (region%peulInput%ptypes(ipt)%negReport == peul_neg_report_none &
        .AND. .NOT.clip) cycle

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

          IF (cv(ipt,ijkc0) < minc) THEN
            minc = cv(ipt,ijkc0)
            iloc1 = i
            jloc1 = j
            kloc1 = k
          ENDIF ! minc

          IF (cv(ipt,ijkc0) > maxc) THEN
            maxc = cv(ipt,ijkc0)
            iloc2 = i
            jloc2 = j
            kloc2 = k
          ENDIF ! maxc

          IF (cv(ipt,ijkc0) < 0._rfreal) THEN
            ncnt1 = ncnt1 + 1
            IF (clip) cv(ipt,ijkc0) = 0._rfreal
          ENDIF ! cv(ipt,ijkC0)

          IF (cv(ipt,ijkc0) > maxconc) THEN
            ncnt2 = ncnt2 + 1
            IF (clip) cv(ipt,ijkc0) = maxconc
          ENDIF ! cv(ipt,ijkC0)

        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k

    IF (ncnt1 > 0 .OR. ncnt2 > 0) THEN

      SELECT CASE (region%peulInput%ptypes(ipt)%negReport)

      CASE (peul_neg_report_used)

        IF (firsttype) THEN

          WRITE(stdout,'(A,1X,A)') solver_name, &
                'Negative smoke concentration detected!'

          IF ( global%flowType == flow_unsteady ) THEN
            WRITE(stdout,'(A,3X,A,1X,ES12.5,A,I1)') solver_name, &
              'Current time and stage:',global%currentTime,', ',region%irkStep
          ELSE
            WRITE(stdout,'(A,3X,A,1X,I6)') solver_name, &
              'Current iteration number:', global%currentIter
          ENDIF ! flowType

          WRITE(stdout,'(A,2X,12(1X,A))') solver_name,    &
                                               ' Region', &
                                                 '  ipt', &
                                               ' # negs', &
                                       '    Least value', &
                                                 '    i', &
                                                 '    j', &
                                                 '    k', &
                                               ' # bigs', &
                                       ' Greatest value', &
                                                 '    i', &
                                                 '    j', &
                                                 '    k'
          firsttype = .false.

        ENDIF ! firstType

        WRITE(stdout,'(A,I10,I6,I8,ES16.6,3I6,I8,ES16.6,3I6)') &
          solver_name,region%iRegionGlobal,ipt,ncnt1,minc,iloc1,jloc1,kloc1, &
          ncnt2,maxc,iloc2,jloc2,kloc2

      END SELECT ! negReport

    ENDIF ! ncnt1 > 0 .OR. ncnt2 > 0

  ENDDO ! ipt

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

  CALL deregisterfunction( global )

END SUBROUTINE peul_enforcepositivity

!******************************************************************************
!
! RCS Revision history:
!
! $Log: PEUL_EnforcePositivity.F90,v $
! Revision 1.3  2008/12/06 08:44:39  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:51  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/12/01 21:09:33  haselbac
! Initial revision after changing case
!
! Revision 1.3  2004/05/03 15:09:05  jferry
! added ability to clip smoke values greater than maximal space packing
!
! Revision 1.2  2004/03/05 22:09:04  jferry
! created global variables for peul, plag, and inrt use
!
! Revision 1.1  2004/03/02 21:44:52  jferry
! Added clipping options
!
!******************************************************************************