RADI_ExtinctionCoef.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: compute diffusion approximation (Rosseland) radiative fluxes by
!          average of variables.
!
! Description: none.
!
! Input: region = data of current region.
!
! Output: region%levels%mixt%rhs = radiative fluxes added to the residual.
!
! Notes: grad(T) is averaged to cell center in wvInt for computation of
!        radiation intensity later.
!
!******************************************************************************
!
! $Id: RADI_ExtinctionCoef.F90,v 1.3 2008/12/06 08:44:37 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE radi_extinctioncoef( region )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getdimensdummy, &
                            rflo_getcelloffset
  USE moderror
  USE modparameters
  USE radi_modparameters
#ifdef PLAG
  USE plag_modparameters
#endif
#ifdef PEUL
  USE peul_modparameters
#endif
  IMPLICIT NONE

#include "Indexing.h"

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

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

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

  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: idcbeg, idcend, jdcbeg, jdcend, kdcbeg, kdcend
  INTEGER :: ilev, icoff, ijcoff, ijkc0, ibc, iec
  INTEGER :: n1, n2, n3, n4, n5, n6, naddgas, naddplag, naddpeul

  REAL(RFREAL) :: rati, maxvfrac, rtresh
  REAL(RFREAL), POINTER :: coef(:,:), optc(:,:), cellvol(:), gofact(:)
  REAL(RFREAL), POINTER :: vfrac(:,:), pdiam(:,:)
  REAL(RFREAL), ALLOCATABLE :: volfrac(:), exteff(:), rdiam(:)
#ifdef PLAG
  INTEGER :: ipcls, npcls, ncont, npclcel
  REAL(RFREAL) :: diasumcel, volsumcel, volsumpcl
  INTEGER,      POINTER :: paiv(:,:), pdvplagvolu(:)
  REAL(RFREAL), POINTER :: pdv(:,:)
#endif
#ifdef PEUL
  INTEGER :: iptype
  REAL(RFREAL) :: ptdiam, pteffdens, smokdens
  REAL(RFREAL), POINTER :: peulcv(:,:)
#endif

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

  CALL registerfunction( region%global,'RADI_ExtinctionCoef',&
  'RADI_ExtinctionCoef.F90' )

! get dimensions, pointers and constants --------------------------------------

  ilev = region%currLevel

  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getdimensdummy( region,ilev,idcbeg,idcend, &
                            jdcbeg,jdcend,kdcbeg,kdcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
  ibc = indijk(idcbeg,jdcbeg,kdcbeg,icoff,ijcoff)
  iec = indijk(idcend,jdcend,kdcend,icoff,ijcoff)

  coef   => region%levels(ilev)%radi%radCoef
  gofact => region%levels(ilev)%radi%goFact
  optc   => region%radiInput%optConst

  rati = 1.5_rfreal

! check phase indexing

  IF ((radi_phase_dispart - radi_phase_gas /= 1 ) .OR. &
      (radi_phase_conpart - radi_phase_dispart /= 1 )) THEN
    CALL errorstop( global,err_radi_fixparam,__line__, &
         'radiation phase indexing is not consistent' )
  ENDIF

! clear up radiation coefficients before accumulating optical constants

  coef = 0._rfreal

! extinction coefficient from contribution of AlOxide smoke and Al-dropplets

  IF (region%radiInput%media == radi_media_artif) THEN

    ALLOCATE( volfrac(nphase), exteff(nphase), rdiam(nphase) )

    DO n = 1,nphase
      volfrac(n) = optc(phase_prop_v,n)
      exteff(n)  = optc(phase_prop_q,n)
      rdiam(n)   = 1._rfreal/optc(phase_prop_d,n)
    ENDDO

    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkc0 = indijk(i ,j ,k ,icoff,ijcoff)

          DO n = 1,nphase
            coef(ijkc0,radi_coeff_extinct)= coef(ijkc0,radi_coeff_extinct)+ &
                                            volfrac(n)*exteff(n)*rdiam(n)
          ENDDO

          coef(ijkc0,radi_coeff_extinct) = rati*coef(ijkc0,radi_coeff_extinct)
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k

    DEALLOCATE( volfrac, exteff, rdiam )

  ELSEIF (region%radiInput%media == radi_media_real) THEN

    cellvol  => region%levels(ilev)%grid%vol
    naddgas  = 0
    naddplag = 0
    naddpeul = 0

#ifdef PLAG
    IF (global%plagUsed) THEN
      naddplag = 0
    ENDIF
#endif
#ifdef PEUL
    IF (global%peulUsed) THEN
      naddpeul = region%peulInput%nPtypes - 1
    ENDIF
#endif

    n1 = radi_phase_gas
    n2 = radi_phase_gas     + naddgas
    n3 = radi_phase_dispart + naddgas
    n4 = radi_phase_dispart + naddgas + naddplag
    n5 = radi_phase_conpart + naddgas + naddplag
    n6 = radi_phase_conpart + naddgas + naddplag + naddpeul

    ALLOCATE( vfrac(n6,ibc:iec) )
    ALLOCATE( pdiam(n6,ibc:iec) )

! - initialiaze volume fractions and particles diameters for safety
    vfrac = 0._rfreal
    pdiam = 1._rfreal
    vfrac(n1:n2,:) = optc(phase_prop_v,radi_phase_gas)
    pdiam(n1:n2,:) = optc(phase_prop_d,radi_phase_gas)

! - disPart volume fraction and diameter
#ifdef PLAG
    IF (global%plagUsed) THEN

      paiv        => region%levels(ilev)%plag%aiv
      pdv         => region%levels(ilev)%plag%dv
      pdvplagvolu => region%levels(ilev)%plag%dvPlagVolu

      npcls    = region%levels(ilev)%plag%nPcls
      ncont    = region%plagInput%nCont

! --- active for non-zero number of particles in region
      IF ( npcls > 0 ) THEN

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

              diasumcel = 0._rfreal
              volsumcel = 0._rfreal
              npclcel   = 0

              DO ipcls = 1, npcls

! ------------- only particles in ijkC0 cell
                IF (paiv(aiv_plag_icells,ipcls) == ijkc0) THEN

! --------------- particle volume
                  volsumpcl = sum(pdv(pdvplagvolu(:),ipcls))

! --------------- volume of all pcls in this cell
                  volsumcel = volsumcel + volsumpcl

! --------------- sum all diameters and number of particles in this cell
                  diasumcel = diasumcel + pdv(dv_plag_diam,ipcls)
                  npclcel   = npclcel + 1
                ENDIF

              ENDDO ! iPcls

              vfrac(n3:n4,ijkc0) = volsumcel/cellvol(ijkc0)
              pdiam(n3:n4,ijkc0) = diasumcel/REAL(npclcel,kind=rfreal)
            ENDDO ! i
          ENDDO   ! j
        ENDDO     ! k

! ----- copy vFrac and pDiam to dummy cells
        CALL copy2dumcells( n3,n4,vfrac )
        CALL copy2dumcells( n3,n4,pdiam )

      ENDIF ! nPcls >0
    ENDIF   ! plagUsed
#endif

! - conPart volume fraction and diameter
#ifdef PEUL
    IF (global%peulUsed) THEN

      peulcv => region%levels(ilev)%peul%cv

! --- there are nPtypes smoke particles, each treated as different phase
      DO iptype = 1,region%peulInput%nPtypes

        ptdiam    = region%peulInput%ptypes(iptype)%diam
        pteffdens = region%peulInput%ptypes(iptype)%denseff

        DO k=kdcbeg,kdcend
          DO j=jdcbeg,jdcend
            DO i=idcbeg,idcend
              ijkc0 = indijk(i ,j ,k ,icoff,ijcoff)

              smokdens = max( peulcv(cv_peul_dens+iptype-1,ijkc0), 0._rfreal )
              vfrac(n5+iptype-1,ijkc0) = smokdens/pteffdens
              pdiam(n5+iptype-1,ijkc0) = ptdiam

            ENDDO ! i
          ENDDO   ! j
        ENDDO     ! k
      ENDDO       ! iPtype
    ENDIF         ! peulUsed
#endif

    ALLOCATE( volfrac(n6) )
    ALLOCATE( exteff(n6) )
    volfrac = 0._rfreal
    exteff  = 0._rfreal

    exteff(n1:n2)  = optc(phase_prop_q,radi_phase_gas)
    exteff(n3:n4)  = optc(phase_prop_q,radi_phase_dispart)
    exteff(n5:n6)  = optc(phase_prop_q,radi_phase_conpart)

    maxvfrac = 0._rfreal

    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkc0 = indijk(i ,j ,k ,icoff,ijcoff)

          DO n = n1, n2
            volfrac(n) = min( optc(phase_prop_v,radi_phase_gas), &
                              vfrac(n,ijkc0) )
          ENDDO
          DO n = n3, n4
            volfrac(n) = min( optc(phase_prop_v,radi_phase_dispart), &
                              vfrac(n,ijkc0) )
          ENDDO
          DO n = n5, n6
            volfrac(n) = min( optc(phase_prop_v,radi_phase_conpart), &
                              vfrac(n,ijkc0) )
          ENDDO

          maxvfrac = max( maxvfrac, sum( volfrac(:) ) )

          DO n = n1,n6
            coef(ijkc0,radi_coeff_extinct)= coef(ijkc0,radi_coeff_extinct)+ &
                                            volfrac(n)*exteff(n)/pdiam(n,ijkc0)
          ENDDO

          coef(ijkc0,radi_coeff_extinct) = rati*coef(ijkc0,radi_coeff_extinct)
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k

    DEALLOCATE( vfrac, pdiam, volfrac, exteff )
    IF (maxvfrac > 1._rfreal) goto 10

  ENDIF         ! media real

! compute go-ahead factor based on EC treshold to obey diffusion approximation

  rtresh = 1._rfreal/radi_real_ecmin

  DO k=kdcbeg,kdcend
    DO j=jdcbeg,jdcend
      DO i=idcbeg,idcend
        ijkc0 = indijk(i ,j ,k ,icoff,ijcoff)
        gofact(ijkc0) = aint( coef(ijkc0,radi_coeff_extinct)*rtresh )
        gofact(ijkc0) = min( gofact(ijkc0), 1._rfreal )
      ENDDO   ! i
    ENDDO     ! j
  ENDDO       ! k

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

  goto 999

! error handling

10  CONTINUE
  CALL errorstop( global,err_radi_mulphase,__line__, &
       'total volume fractions of particles and smoke > 1.' )

999 CONTINUE

  CALL deregisterfunction( region%global )

! =============================================================================
!   Copy given vector variables to dummy cells
! =============================================================================

CONTAINS

  SUBROUTINE copy2dumcells( m1,m2,var )

! ... parameters
    INTEGER :: m1, m2, ibegc, iendc
    REAL(RFREAL), POINTER :: var(:,:)

! ... local variables
    INTEGER :: l, ijkc, ijkci

! - I direction

    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend

        DO i=idcbeg,ipcbeg-1
          ijkc  = indijk(i     ,j     ,k     ,icoff,ijcoff)
          ijkci = indijk(1     ,j     ,k     ,icoff,ijcoff)

          DO l = m1, m2
            var(l,ijkc) = var(l,ijkci)
          ENDDO
        ENDDO
        DO i=ipcend+1,idcend
          ijkc  = indijk(i      ,j     ,k     ,icoff,ijcoff)
          ijkci = indijk(ipcend ,j     ,k     ,icoff,ijcoff)

          DO l = m1, m2
            var(l,ijkc) = var(l,ijkci)
          ENDDO
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k

! - J direction

    DO k=kdcbeg,kdcend
      DO i=idcbeg,idcend

        DO j=jdcbeg,jpcbeg-1
          ijkc  = indijk(i     ,j     ,k     ,icoff,ijcoff)
          ijkci = indijk(i     ,1     ,k     ,icoff,ijcoff)

          DO l = m1, m2
            var(l,ijkc) = var(l,ijkci)
          ENDDO
        ENDDO

        DO j=jpcend+1,jdcend
          ijkc  = indijk(i     ,j      ,k     ,icoff,ijcoff)
          ijkci = indijk(i     ,jpcend ,k     ,icoff,ijcoff)

          DO l = m1, m2
            var(l,ijkc) = var(l,ijkci)
          ENDDO
        ENDDO   ! j
      ENDDO     ! i
    ENDDO       ! k

! - K direction

    DO j=jdcbeg,jdcend
      DO i=idcbeg,idcend

        DO k=kdcbeg,kpcbeg-1
          ijkc  = indijk(i     ,j     ,k     ,icoff,ijcoff)
          ijkci = indijk(i     ,j     ,1     ,icoff,ijcoff)

          DO l = m1, m2
            var(l,ijkc) = var(l,ijkci)
          ENDDO
        ENDDO

        DO k=kpcend+1,kdcend
          ijkc  = indijk(i     ,j     ,k      ,icoff,ijcoff)
          ijkci = indijk(i     ,j     ,kpcend ,icoff,ijcoff)

          DO l = m1, m2
            var(l,ijkc) = var(l,ijkci)
          ENDDO
        ENDDO   ! k
      ENDDO     ! i
    ENDDO       ! j

  END SUBROUTINE copy2dumcells

END SUBROUTINE radi_extinctioncoef

!******************************************************************************
!
! RCS Revision history:
!
! $Log: RADI_ExtinctionCoef.F90,v $
! Revision 1.3  2008/12/06 08:44:37  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:49  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/09/22 02:35:49  wasistho
! changed file nomenclature from lower to upper case
!
! Revision 1.3  2004/03/05 22:09:04  jferry
! created global variables for peul, plag, and inrt use
!
! Revision 1.2  2003/07/30 22:23:02  wasistho
! enter part and smoke data into radiation
!
! Revision 1.1  2003/07/17 01:16:59  wasistho
! initial activation rocrad
!
!
!******************************************************************************