RADI_FluxLimiter.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 the radiation flux limiter based on the model selected.
!
! Description: none.
!
! Input: region = data of current region.
!
! Output: region%levels%radi%fluxLim = flux limiter in diff. based radiation
!
! Notes: none. 
!
!******************************************************************************
!
! $Id: RADI_FluxLimiter.F90,v 1.3 2008/12/06 08:44:37 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE radi_fluxlimiter( region )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modinterfaces, ONLY : rflo_getdimensdummy, rflo_getcelloffset, &
                            rflo_getnodeoffset
  USE moderror
  USE modparameters
  USE radi_modparameters
  IMPLICIT NONE

#include "Indexing.h"

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

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

! ... local variables
  INTEGER :: idcbeg, idcend, jdcbeg, jdcend, kdcbeg, kdcend
  INTEGER :: ilev, icoff, ijcoff, inoff, ijnoff
  INTEGER :: ijkc, ijkn, ijkni, ijknj, ijknk, errfl

  REAL(RFREAL) :: one3rd, one6th, rparam, grademod, coefc, radenerg
  REAL(RFREAL), POINTER :: dv(:,:), rcv(:,:), rdv(:,:), coef(:,:), fluxlim(:) 
  REAL(RFREAL), POINTER :: gradi(:,:), gradj(:,:), gradk(:,:)
  REAL(RFREAL), POINTER :: rgradi(:,:), rgradj(:,:), rgradk(:,:)
  REAL(RFREAL), ALLOCATABLE :: gradc(:,:,:,:)

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

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

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

  one3rd = 1._rfreal/3._rfreal
  one6th = 1._rfreal/6._rfreal
  ilev   = region%currLevel

  CALL rflo_getdimensdummy( region,ilev,idcbeg,idcend, &
                            jdcbeg,jdcend,kdcbeg,kdcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )

  gradi  => region%levels(ilev)%mixt%gradi
  gradj  => region%levels(ilev)%mixt%gradj
  gradk  => region%levels(ilev)%mixt%gradk
  dv     => region%levels(ilev)%mixt%dv
  rdv    => region%levels(ilev)%radi%dv
  coef   => region%levels(ilev)%radi%radCoef

  IF (region%radiInput%radiModel == radi_model_fldsrc .OR. & 
      region%radiInput%radiModel == radi_model_fldtran) THEN
    ALLOCATE( gradc(idcbeg:idcend,jdcbeg:jdcend,kdcbeg:kdcend,3),stat=errfl )
    region%global%error = errfl
    CALL errorstop( region%global,err_allocate,__line__ )
  ENDIF
 
! get effective temperature --------------------------------------------------

  IF ((region%radiInput%radiModel == radi_model_ross) .OR. &
      (region%radiInput%radiModel == radi_model_fldsrc .AND. & 
       region%radiInput%fluxLim == fld_lim_none)) THEN
    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkc = indijk(i  ,j,k,icoff,ijcoff)
          fluxlim(ijkc) = one3rd
        ENDDO   ! i
      ENDDO   ! j
    ENDDO   ! k
  ELSEIF ((region%radiInput%radiModel == radi_model_fldsrc) .AND. &
          (region%radiInput%fluxLim == fld_lim_lp)) THEN

    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkn  = indijk(i  ,j  ,k  ,inoff,ijnoff)
          ijkni = indijk(i+1,j  ,k  ,inoff,ijnoff)
          ijknj = indijk(i  ,j+1,k  ,inoff,ijnoff)
          ijknk = indijk(i  ,j  ,k+1,inoff,ijnoff)
          gradc(i,j,k,1) = one6th* &
                         (gradi(gr_mixt_tx,ijkn)+gradi(gr_mixt_tx,ijkni) + &
                          gradj(gr_mixt_tx,ijkn)+gradj(gr_mixt_tx,ijknj) + &
                          gradk(gr_mixt_tx,ijkn)+gradk(gr_mixt_tx,ijknk))
          gradc(i,j,k,2) = one6th* &
                         (gradi(gr_mixt_ty,ijkn)+gradi(gr_mixt_ty,ijkni) + &
                          gradj(gr_mixt_ty,ijkn)+gradj(gr_mixt_ty,ijknj) + &
                          gradk(gr_mixt_ty,ijkn)+gradk(gr_mixt_ty,ijknk))
          gradc(i,j,k,3) = one6th* &
                         (gradi(gr_mixt_tz,ijkn)+gradi(gr_mixt_tz,ijkni) + &
                          gradj(gr_mixt_tz,ijkn)+gradj(gr_mixt_tz,ijknj) + &
                          gradk(gr_mixt_tz,ijkn)+gradk(gr_mixt_tz,ijknk))
        ENDDO   ! i
      ENDDO   ! j
    ENDDO   ! k
    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkc = indijk(i  ,j,k,icoff,ijcoff)
          grademod = sqrt( gradc(i,j,k,1)*gradc(i,j,k,1) + &
                           gradc(i,j,k,2)*gradc(i,j,k,2) + &
                           gradc(i,j,k,3)*gradc(i,j,k,3) )
          coefc    = coef(ijkc,radi_coeff_extinct)
          radenerg = abs( dv(dv_mixt_temp,ijkc) )
          rparam   = grademod/(coefc*radenerg)
          fluxlim(ijkc) = flimiterlp(rparam)
        ENDDO   ! i
      ENDDO   ! j
    ENDDO   ! k
  ELSEIF (region%radiInput%radiModel == radi_model_fldtran) THEN

    rcv    => region%levels(ilev)%radi%cv
    rgradi => region%levels(ilev)%radi%gradi
    rgradj => region%levels(ilev)%radi%gradj
    rgradk => region%levels(ilev)%radi%gradk

    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkn  = indijk(i  ,j  ,k  ,inoff,ijnoff)
          ijkni = indijk(i+1,j  ,k  ,inoff,ijnoff)
          ijknj = indijk(i  ,j+1,k  ,inoff,ijnoff)
          ijknk = indijk(i  ,j  ,k+1,inoff,ijnoff)
          gradc(i,j,k,1) = one6th* &
                         (rgradi(gr_radi_ex,ijkn)+rgradi(gr_radi_ex,ijkni) + &
                          rgradj(gr_radi_ex,ijkn)+rgradj(gr_radi_ex,ijknj) + &
                          rgradk(gr_radi_ex,ijkn)+rgradk(gr_radi_ex,ijknk))
          gradc(i,j,k,2) = one6th* &
                         (rgradi(gr_radi_ey,ijkn)+rgradi(gr_radi_ey,ijkni) + &
                          rgradj(gr_radi_ey,ijkn)+rgradj(gr_radi_ey,ijknj) + &
                          rgradk(gr_radi_ey,ijkn)+rgradk(gr_radi_ey,ijknk))
          gradc(i,j,k,3) = one6th* &
                         (rgradi(gr_radi_ez,ijkn)+rgradi(gr_radi_ez,ijkni) + &
                          rgradj(gr_radi_ez,ijkn)+rgradj(gr_radi_ez,ijknj) + &
                          rgradk(gr_radi_ez,ijkn)+rgradk(gr_radi_ez,ijknk))
        ENDDO   ! i
      ENDDO   ! j
    ENDDO   ! k

    DO k=kdcbeg,kdcend
      DO j=jdcbeg,jdcend
        DO i=idcbeg,idcend
          ijkc = indijk(i  ,j,k,icoff,ijcoff)
          grademod = sqrt( gradc(i,j,k,1)*gradc(i,j,k,1) + &
                           gradc(i,j,k,2)*gradc(i,j,k,2) + &
                           gradc(i,j,k,3)*gradc(i,j,k,3) )
          coefc    = coef(ijkc,radi_coeff_extinct)
          radenerg = rcv(cv_radi_ener,ijkc)
          rparam   = grademod/(coefc*radenerg)
          fluxlim(ijkc) = flimiterlp(rparam)
        ENDDO   ! i
      ENDDO   ! j
    ENDDO   ! k
  ENDIF

! deallocate arrays

  IF (ALLOCATED( gradc )) THEN
    DEALLOCATE( gradc, stat=errfl )
    region%global%error = errfl
    CALL errorstop( region%global,err_deallocate,__line__ )
  ENDIF

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

  CALL deregisterfunction( region%global )

! #############################################################################

  CONTAINS

    REAL(RFREAL) FUNCTION flimiterlp( r )

      REAL(RFREAL) :: r

      flimiterlp = (2._rfreal + r)/(6._rfreal + 3._rfreal*r + r*r)

    END FUNCTION flimiterlp

END SUBROUTINE radi_fluxlimiter

!******************************************************************************
!
! RCS Revision history:
!
! $Log: RADI_FluxLimiter.F90,v $
! Revision 1.3  2008/12/06 08:44:37  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:50  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.2  2004/09/22 01:32:19  wasistho
! switch LFD to FLD for flux limited diffusion
!
! Revision 1.1  2004/09/18 18:02:07  wasistho
! initial import for Limited Flux Diffusion radiation
!
!
!
!******************************************************************************