PERI_coCprSlowTermsFlo.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 additional terms to NS due to CPR formulations
!
! Description: all variables in slow terms are averaged over i,k plane.
!
! Input: region = data of current region.
!
! Output: region%levels%mixt%rhs = CPR slow terms added to the residual.
!
! Notes: This routine contents MPI global SUM which does not account for
!        summation over regions in same processor. Therefore nProcAlloc should
!        be equal to nRegions.
!
!******************************************************************************
!
! $Id: PERI_coCprSlowTermsFlo.F90,v 1.3 2008/12/06 08:44:37 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE peri_cocprslowterms( region )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
                            mixtperf_r_m, mixtperf_g_cpr
  USE moderror
  USE modmpi
  USE modparameters
  USE peri_modparameters
  IMPLICIT NONE

#include "Indexing.h"

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

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

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

  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ilev, icoff, ijcoff, ijkc0, nvar, ipc, jpc, kpc
  INTEGER :: indcp, indmol

  REAL(RFREAL), POINTER :: cv(:,:), dv(:,:), gv(:,:), rhs(:,:), vol(:)
  REAL(RFREAL), POINTER :: cprvar(:,:), varsend(:,:), varrecv(:,:)
  REAL(RFREAL), ALLOCATABLE :: rgas(:), cpgas(:)
  REAL(RFREAL) :: voleps, aversize, sndaversize, rcvaversize, ravgrho
  REAL(RFREAL) :: denumer, denomin, delta, rdelta, gamma, gammin, gagmin
  REAL(RFREAL) :: rho, ruc, ren, uve, vve, wve, tmp, prs, rprim, uprim, vprim 
  REAL(RFREAL) :: wprim, tprim, pprim, ugrad, vgrad, wgrad, tgrad, pgrad
  REAL(RFREAL) :: refmeanpgrad, fcont, fmomx, fmomy, fmomz, frhoe

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

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

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

! get dimensions, pointers and parameters -------------------------------------

  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
  dv     => region%levels(ilev)%mixt%dv
  gv     => region%levels(ilev)%mixt%gv
  rhs    => region%levels(ilev)%mixt%rhs
  vol    => region%levels(ilev)%grid%vol
  cprvar => region%levels(ilev)%peri%cprVar
  varsend=> region%levels(ilev)%peri%varSend
  varrecv=> region%levels(ilev)%peri%varRecv

  ipc      = ipcend-ipcbeg+1
  jpc      = jpcend-jpcbeg+1
  kpc      = kpcend-kpcbeg+1
  aversize = dble(ipc*kpc)

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

  refmeanpgrad = region%periInput%meanPgrad
  delta        = region%global%refLength
  rdelta       = 1._rfreal/delta

! allocate temporary space

  ALLOCATE( rgas(jpcbeg:jpcend), cpgas(jpcbeg:jpcend) )

! x-z plane averaged variables

  DO j=jpcbeg,jpcend
    cprvar(1:nvar,j) = 0._rfreal
    rgas(j)          = 0._rfreal
    cpgas(j)         = 0._rfreal
    DO k=kpcbeg,kpcend
      DO i=ipcbeg,ipcend
        ijkc0 = indijk(i ,j ,k ,icoff,ijcoff)

        cprvar(cpr_rho,j) = cprvar(cpr_rho,j) + cv(cv_mixt_dens,ijkc0)
        cprvar(cpr_ruc,j) = cprvar(cpr_ruc,j) + cv(cv_mixt_xmom,ijkc0)
        cprvar(cpr_rvc,j) = cprvar(cpr_rvc,j) + cv(cv_mixt_ymom,ijkc0)
        cprvar(cpr_uve,j) = cprvar(cpr_uve,j) + dv(dv_mixt_uvel,ijkc0)
        cprvar(cpr_vve,j) = cprvar(cpr_vve,j) + dv(dv_mixt_vvel,ijkc0)
        cprvar(cpr_tmp,j) = cprvar(cpr_tmp,j) + dv(dv_mixt_temp,ijkc0)
        cprvar(cpr_prs,j) = cprvar(cpr_prs,j) + dv(dv_mixt_pres,ijkc0)

        rgas(j) = rgas(j) + mixtperf_r_m( gv(gv_mixt_mol,ijkc0*indmol) )
        cpgas(j)= cpgas(j)+ gv(gv_mixt_cp,ijkc0*indcp)
      ENDDO
    ENDDO
  ENDDO

#ifdef MPI
  DO j = jpcbeg, jpcend
    DO l = 1,nvar
      varsend(j-jpcbeg+1,l) = cprvar(l,j)
    END DO
    varsend(j-jpcbeg+1,nvar+1) = rgas(j)
    varsend(j-jpcbeg+1,nvar+2) = cpgas(j)
  END DO
  CALL mpi_barrier( region%global%mpiComm, region%global%mpierr )
  DO l = 1, nvar+gas_nvar
    CALL mpi_allreduce( varsend(1,l),varrecv(1,l),jpc,mpi_double_precision, &
         mpi_sum,region%global%mpiComm, region%global%mpierr )
  END DO
  DO j = jpcbeg, jpcend
    DO l = 1,nvar
      cprvar(l,j) = varrecv(j-jpcbeg+1,l)
    END DO
    rgas(j)  = varrecv(j-jpcbeg+1,nvar+1)
    cpgas(j) = varrecv(j-jpcbeg+1,nvar+2)
  END DO

  sndaversize = aversize
  CALL mpi_allreduce( sndaversize,rcvaversize,1, mpi_double_precision, &
       mpi_sum,region%global%mpiComm,region%global%mpierr )
  aversize = rcvaversize
#endif

  DO j = jpcbeg,jpcend
    DO l = 1,nvar
      cprvar(l,j) = cprvar(l,j)/aversize
    ENDDO
    rgas(j) = rgas(j)/aversize
    cpgas(j)= cpgas(j)/aversize

    gamma   = mixtperf_g_cpr( cpgas(j),rgas(j) )
    gagmin  = gamma/(gamma - 1_rfreal)

    ravgrho = 1._rfreal/cprvar(cpr_rho,j)  
    denumer = gagmin*refmeanpgrad + (cprvar(cpr_ruc,j)**2 + &
                                     cprvar(cpr_rvc,j)**2)*rdelta*ravgrho
    denomin = gagmin*cprvar(cpr_prs,j)*ravgrho + &
              (cprvar(cpr_ruc,j)**2+cprvar(cpr_rvc,j)**2)*ravgrho**2
    cprvar(cpr_dor,j) = denumer/denomin
    cprvar(cpr_dou,j) = (cprvar(cpr_ruc,j)*rdelta-cprvar(cpr_uve,j)* &
                        cprvar(cpr_dor,j))*ravgrho
    cprvar(cpr_dov,j) = -cprvar(cpr_vve,j)*ravgrho*cprvar(cpr_dor,j)
    cprvar(cpr_dot,j) = (refmeanpgrad/rgas(j)-cprvar(cpr_tmp,j)* &
                        cprvar(cpr_dor,j))*ravgrho
  ENDDO

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

        rho = cv(cv_mixt_dens,ijkc0)
        ruc = cv(cv_mixt_xmom,ijkc0)
        ren = cv(cv_mixt_ener,ijkc0)
        uve = dv(dv_mixt_uvel,ijkc0)
        vve = dv(dv_mixt_vvel,ijkc0)
        wve = dv(dv_mixt_wvel,ijkc0)
        tmp = dv(dv_mixt_temp,ijkc0)
        prs = dv(dv_mixt_pres,ijkc0)

        rprim = rho - cprvar(cpr_rho,j)
        uprim = uve - cprvar(cpr_uve,j)
        vprim = vve - cprvar(cpr_vve,j)
        wprim = wve
        tprim = tmp - cprvar(cpr_tmp,j)
        pprim = prs - cprvar(cpr_prs,j)

        ugrad = uprim*rdelta + cprvar(cpr_dou,j)
        vgrad = vprim*rdelta + cprvar(cpr_dov,j)  
        wgrad = wprim*rdelta
        tgrad = tprim*rdelta + cprvar(cpr_dot,j)  
        pgrad = pprim*rdelta + refmeanpgrad 

        fcont = (cprvar(cpr_ruc,j)+cprvar(cpr_rho,j)*uprim+cprvar(cpr_uve,j)* &
                rprim)*rdelta + uprim*cprvar(cpr_dor,j) + rprim*cprvar(cpr_dou,j)
        fmomx = uve*fcont + ruc*ugrad + pgrad
        fmomy = vve*fcont + ruc*vgrad
        fmomz = wve*fcont + ruc*wgrad
        gamma = mixtperf_g_cpr( cpgas(j),rgas(j) )
        gammin= gamma-1._rfreal
        frhoe = ren/rho*fcont + ruc*(uve*ugrad+vve*vgrad+wve*wgrad+rgas(j)/ &
                gammin*tgrad) + prs*ugrad +uve*pgrad

        voleps= vol(ijkc0)*region%periInput%cprEpsilon
        rhs(cv_mixt_dens,ijkc0) = rhs(cv_mixt_dens,ijkc0) + voleps*fcont
        rhs(cv_mixt_xmom,ijkc0) = rhs(cv_mixt_xmom,ijkc0) + voleps*fmomx
        rhs(cv_mixt_ymom,ijkc0) = rhs(cv_mixt_ymom,ijkc0) + voleps*fmomy
        rhs(cv_mixt_zmom,ijkc0) = rhs(cv_mixt_zmom,ijkc0) + voleps*fmomz
        rhs(cv_mixt_ener,ijkc0) = rhs(cv_mixt_ener,ijkc0) + voleps*frhoe
!=======================================
!        fCont = cprVar(CPR_RUC,j)*rdelta
!        fMomX = cprVar(CPR_UVE,j)*fCont 
!        fMomY = cprVar(CPR_VVE,j)*fCont
!        fMomZ = 0._RFREAL
!        fRhoE = ren/rho*fCont + ruc*(uve*uGrad+vve*vGrad+wve*wGrad+rgas(j)/ &
!                gammin*tGrad) + prs*uGrad +uve*pGrad
!        rhs(CV_MIXT_DENS,ijkC0) = rhs(CV_MIXT_DENS,ijkC0) + volEps*fCont
!        rhs(CV_MIXT_XMOM,ijkC0) = rhs(CV_MIXT_XMOM,ijkC0) + volEps* &
!                                  (fMomX + ruc*uve*rdelta + refMeanPgrad)
!        rhs(CV_MIXT_YMOM,ijkC0) = rhs(CV_MIXT_YMOM,ijkC0) + volEps*ruc*vve*rdelta
!        rhs(CV_MIXT_ZMOM,ijkC0) = rhs(CV_MIXT_ZMOM,ijkC0) + volEps*ruc*wve*rdelta
!        rhs(CV_MIXT_ENER,ijkC0) = rhs(CV_MIXT_ENER,ijkC0) + volEps* &
!                                  (ren/rho+prs)*uve*rdelta
!========================================
      ENDDO
    ENDDO
  ENDDO

! deallocate temporary workspace

  DEALLOCATE( rgas, cpgas )

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

  CALL deregisterfunction( region%global )

END SUBROUTINE peri_cocprslowterms

!******************************************************************************
!
! RCS Revision history:
!
! $Log: PERI_coCprSlowTermsFlo.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/06/08 23:56:56  wasistho
! changed nomenclature
!
! Revision 1.2  2003/05/15 02:57:05  jblazek
! Inlined index function.
!
! Revision 1.1.1.1  2003/03/29 03:36:30  wasistho
! install ROCPERI
!
!
!******************************************************************************