PERI_coPgradUpdateFlo.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: update pressure gradient based on mass flux balance
!
! Description: when computed mass flux is lower than reference mass flux
!              pressure gradient is adjusted in such away that the bulk flow 
!              accelerate, in other words mass flux is conserved
!
! Input: region = data of current region.
!
! Output: region%levels%mixt%rhs = updated pressure gradient added to residual.
!
! Notes: If MPI, nProcAlloc should be equal to nRegions, otherwise the wall
!        normal extent should be covered in one region.
!
!******************************************************************************
!
! $Id: PERI_coPgradUpdateFlo.F90,v 1.8 2008/12/06 08:44:37 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE peri_copgradupdate( region )

  USE moddatatypes
  USE modglobal, ONLY     : t_global
  USE moddatastruct, ONLY : t_region
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
                            rflo_getnodeoffset
  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, n, ir

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

  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ilev, icoff, ijcoff, inoff, ijnoff, ijkc0, ijkc1
  INTEGER :: ijkn, ijknp1, ipc, kpc, istage

  REAL(RFREAL)  :: rhoc, rhouc, rhovc, rhowc, rhoec, pc
  REAL(RFREAL)  :: refmassflux, rescaler, alpha, massflux
  REAL(RFREAL)  :: cfl, dt, rucm, delta, deltamassflux, aversize
  REAL(RFREAL)  :: sndmassflux, rcvmassflux, sndaversize,  rcvaversize
  REAL(RFREAL)  :: snddt, rcvdt, avgdt, dtfactor, refmeanpgrad, reftauwall
  REAL(RFREAL)  :: tauwall, avgtauwall, sndtauwall, rcvtauwall
  REAL(RFREAL)  :: rn, sndrn, rcvrn
  REAL(RFREAL), POINTER :: xyz(:,:), cv(:,:), dv(:,:), tv(:,:), rdt(:)

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

  rcsidentstring = '$RCSfile: PERI_coPgradUpdateFlo.F90,v $ $Revision: 1.8 $'

  global => region%global
  CALL registerfunction( global,'PERI_CoPgradUpdate',&
  'PERI_coPgradUpdateFlo.F90' )

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

  ilev = region%currLevel
  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
 
  xyz => region%levels(ilev)%grid%xyz
  cv  => region%levels(ilev)%mixt%cv
  dv  => region%levels(ilev)%mixt%dv
  tv  => region%levels(ilev)%mixt%tv
  rdt => region%levels(ilev)%dt

  istage      = mod(region%irkStep,global%nrkSteps) + 1
  alpha       = region%mixtInput%ark(istage)
  refmassflux = region%periInput%bulkmFlux

  IF (region%periInput%flowKind == peri_flow_cpr) THEN
    rescaler = region%periInput%cprEpsilon
    IF (global%flowType == flow_unsteady .AND. &
        global%solverType == solv_explicit ) alpha = 1._rfreal
  ELSEIF (region%periInput%flowKind == peri_flow_channel) THEN
    rescaler = 1._rfreal
  ENDIF

  IF (global%flowType==flow_steady) THEN
    cfl   = region%mixtInput%cfl
    avgdt = 0._rfreal
    IF (global%solverType == solv_implicit ) THEN ! Dual Tst
      dtfactor = 10._rfreal
    ELSE
!      dtFactor = 0.9_RFREAL
      dtfactor = 10._rfreal
    ENDIF
    n = 0
    DO k=kpcbeg,kpcend
      DO j=jpcbeg,jpcend
        DO i=ipcbeg,ipcend
          n = n+1
          ijkc0 = indijk(i ,j  ,k ,icoff,ijcoff)
          avgdt = avgdt + rdt(ijkc0)
        ENDDO
      ENDDO
    ENDDO
    rn    = REAL(n)
#ifdef MPI
    snddt = avgdt
    CALL mpi_allreduce( snddt,rcvdt,1,mpi_double_precision, mpi_sum, &
                        global%mpiComm, global%mpierr )
    avgdt = rcvdt
    sndrn = rn
    CALL mpi_allreduce( sndrn,rcvrn,1,mpi_double_precision, mpi_sum, &
                        global%mpiComm, global%mpierr )
    rn = rcvrn
#endif
    avgdt = dtfactor*avgdt/rn
    dt    = cfl*avgdt
!    IF (global%solverType == SOLV_IMPLICIT ) dt = global%dtMin ! Dual Tst
  ELSE
    dt    = global%dtMin
  ENDIF

! check if pressure gradient constant in all regions

  refmeanpgrad = region%periInput%meanPgrad

#ifdef MPI
  IF (refmeanpgrad /= global%moduleVar(1) ) THEN
    CALL errorstop( global,err_peri_physparam,__line__, &
         'mean pressure gradient vary between regions' )
  ENDIF
#endif 

! compute mass flux prior rho.u update -------------------------------------

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

  massflux = 0._rfreal

  DO j = jpcbeg,jpcend
    rucm = 0._rfreal
    DO k=kpcbeg,kpcend
      DO i=ipcbeg,ipcend
        ijkc0 = indijk(i ,j  ,k ,icoff,ijcoff)
        ijkn  = indijk(i ,j  ,k ,inoff,ijnoff)
        ijknp1= indijk(i ,j+1,k ,inoff,ijnoff)
        rucm = rucm + cv(cv_mixt_xmom,ijkc0)
      END DO
    END DO
    massflux = massflux + rucm*(xyz(ycoord,ijknp1)-xyz(ycoord,ijkn))
  END DO ! j  

! this part only for channel flow

  avgtauwall = 0._rfreal

  IF (region%periInput%flowKind == peri_flow_channel .AND. &
      region%periInput%split(jcoord) /= off ) THEN

    IF (region%procId == masterproc) THEN
      DO k=kpcbeg,kpcend
        DO i=ipcbeg,ipcend
          ijkc0 = indijk(i ,1  ,k ,icoff,ijcoff)
          ijkc1 = indijk(i ,2  ,k ,icoff,ijcoff)
          ijkn  = indijk(i ,1  ,k ,inoff,ijnoff)
          ijknp1= indijk(i ,2  ,k ,inoff,ijnoff)
          tauwall = 0.5_rfreal*tv(tv_mixt_muel,ijkc0)* &
                    (dv(dv_mixt_uvel,ijkc1)+dv(dv_mixt_uvel,ijkc0))/ &
                    (xyz(ycoord,ijknp1)-xyz(ycoord,ijkn))
          avgtauwall = avgtauwall + tauwall
        ENDDO
      ENDDO
    ENDIF  

    IF (region%procId == global%nProcAlloc-1) THEN
      DO k=kpcbeg,kpcend
        DO i=ipcbeg,ipcend
          ijkc0 = indijk(i ,jpcend   ,k ,icoff,ijcoff)
          ijkc1 = indijk(i ,jpcend-1 ,k ,icoff,ijcoff)
          ijkn  = indijk(i ,jpcend   ,k ,inoff,ijnoff)
          ijknp1= indijk(i ,jpcend+1 ,k ,inoff,ijnoff)
          tauwall = 0.5_rfreal*tv(tv_mixt_muel,ijkc0)* &
                    (dv(dv_mixt_uvel,ijkc1)+dv(dv_mixt_uvel,ijkc0))/ &
                    abs( xyz(ycoord,ijknp1)-xyz(ycoord,ijkn) )
          avgtauwall = avgtauwall + tauwall
        ENDDO
      ENDDO
    ENDIF   ! procId
  ENDIF     ! flowKind=channel

! get average mass flux over all processors

#ifdef MPI
  sndmassflux = massflux
  CALL mpi_allreduce( sndmassflux,rcvmassflux,1,mpi_double_precision,mpi_sum, &
                      global%mpiComm, global%mpierr )
  massflux = rcvmassflux

  IF (region%periInput%flowKind == peri_flow_channel) THEN
    sndtauwall = avgtauwall
    CALL mpi_allreduce( sndtauwall,rcvtauwall,1,mpi_double_precision,mpi_sum, &
                        global%mpiComm, global%mpierr )
    avgtauwall = rcvtauwall
  ENDIF

  IF (region%periInput%split(jcoord) == off) THEN
    sndaversize = aversize
    CALL mpi_allreduce( sndaversize,rcvaversize,1,mpi_double_precision,mpi_sum, &
                        global%mpiComm, global%mpierr )
    aversize = rcvaversize
  ENDIF
#endif
  IF (global%myProcId==masterproc) write(*,*)massflux,aversize

! channel half width needed to update pressure gradient

  delta      = global%refLength
  massflux   = massflux/aversize
  avgtauwall = 0.5_rfreal*avgtauwall/aversize

! update pressure gradient

  deltamassflux = refmassflux - massflux
  IF (global%myProcId==masterproc) WRITE(*,*) region%localNumber, &
    'refMassFlux-massFlux prior updating u',refmassflux,massflux,deltamassflux, &
                                            avgtauwall/delta

  refmeanpgrad  = refmeanpgrad - 0.5_rfreal*deltamassflux/ &
                  (delta*rescaler*alpha*dt)

  IF ((region%periInput%flowKind == peri_flow_channel) .AND. &
      (region%periInput%pgradType == cnl_pgrad_tauwall)) THEN
!    refTauWall   = global%refDensity*region%periInput%cnlUtau**2
!    refMeanPgrad = -refTauWall/global%refLength
    refmeanpgrad = min( refmeanpgrad, -avgtauwall/global%refLength )
  ENDIF

  IF (global%myProcId==masterproc) WRITE(*,*) region%localNumber, &
    ' meanPgrad ',refmeanpgrad,dt,alpha

! store pressure gradient in data structure
  region%periInput%meanPgrad = refmeanpgrad

! update x-momentum due to updated pressure gradient/ mass flux ---------------

  IF (region%periInput%flowKind == peri_flow_cpr) THEN

    massflux = 0._rfreal

    DO j = jpcbeg,jpcend
      rucm = 0._rfreal
      DO k=kpcbeg,kpcend
        DO i=ipcbeg,ipcend
          ijkc0 = indijk(i ,j  ,k ,icoff,ijcoff)
          ijkn  = indijk(i ,j  ,k ,inoff,ijnoff)
          ijknp1= indijk(i ,j+1,k ,inoff,ijnoff)
          cv(cv_mixt_xmom,ijkc0) = cv(cv_mixt_xmom,ijkc0) + 0.5_rfreal* &
                                   deltamassflux/delta
          dv(dv_mixt_uvel,ijkc0) = cv(cv_mixt_xmom,ijkc0)/cv(cv_mixt_dens,ijkc0)
          rucm = rucm + cv(cv_mixt_xmom,ijkc0)
        ENDDO
      ENDDO
      massflux = massflux + rucm*(xyz(ycoord,ijknp1)-xyz(ycoord,ijkn))
    ENDDO ! j  

#ifdef MPI
    CALL mpi_barrier( global%mpiComm, global%mpierr )
    sndmassflux = massflux
    CALL mpi_allreduce( sndmassflux,rcvmassflux,1,mpi_double_precision,mpi_sum, &
                        global%mpiComm, global%mpierr )
    massflux = rcvmassflux
#endif

    massflux = massflux/aversize

    deltamassflux = refmassflux - massflux
    IF (global%myProcId==masterproc) WRITE(*,*) region%localNumber, &
      'refMassFlux-massFlux AFTER updating u',refmassflux,massflux,deltamassflux

  ENDIF ! flowKind=CPR

! save pressure gradient in a global variable for restart purpose

  global%moduleVar(1) = refmeanpgrad

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

  CALL deregisterfunction( global )

END SUBROUTINE peri_copgradupdate

!******************************************************************************
!
! RCS Revision history:
!
! $Log: PERI_coPgradUpdateFlo.F90,v $
! Revision 1.8  2008/12/06 08:44:37  mtcampbe
! Updated license.
!
! Revision 1.7  2008/11/19 22:17:49  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.6  2005/03/08 06:00:33  wasistho
! set dtFactor to 10 for steady computation
!
! Revision 1.5  2004/12/11 03:48:18  wasistho
! resereve dt=global%dtMin for dual tst
!
! Revision 1.4  2004/12/08 19:42:07  wasistho
! used dtFactor
!
! Revision 1.3  2004/12/08 01:47:46  wasistho
! bug fix, MPI_MAX to MPI_SUM
!
! Revision 1.2  2004/12/04 06:49:14  wasistho
! take global%dtMin for dt regardless steady or unsteady
!
! Revision 1.1  2004/06/08 23:56:56  wasistho
! changed nomenclature
!
! Revision 1.15  2004/01/23 22:42:49  wasistho
! set avgDt to 0.9 of the true value for stability of steady state computations
!
! Revision 1.14  2004/01/23 03:11:22  wasistho
! modified dt, added computed tauwall, unst.impl.condition, no cnl massflx correction
!
! Revision 1.13  2003/10/21 20:26:37  wasistho
! global max for dt of unsteady flow
!
! Revision 1.12  2003/10/20 00:32:07  wasistho
! modified istage definition
!
! Revision 1.11  2003/10/18 00:47:50  wasistho
! modified dt for unsteady flow
!
! Revision 1.10  2003/09/18 01:56:26  wasistho
! added ijksplit and pgradType in PERI_PgradUpdate
!
! Revision 1.9  2003/09/12 20:57:45  wasistho
! update pressure gradient only for CPR, fixed for channel
!
! Revision 1.8  2003/09/07 23:32:02  wasistho
! shift ISTAGE 1 stage forward to be consistent
!
! Revision 1.7  2003/05/15 02:57:05  jblazek
! Inlined index function.
!
! Revision 1.6  2003/04/25 23:17:49  wasistho
! update cpr pgrad per time-step for unsteady
!
! Revision 1.5  2003/04/05 04:06:39  wasistho
! shift forward time stage
!
! Revision 1.4  2003/04/04 22:27:19  wasistho
! correct averSize for multi processors
!
! Revision 1.3  2003/04/03 20:57:02  wasistho
! replace regions to region in pgradUpdate
!
! Revision 1.2  2003/04/03 00:30:07  wasistho
! uniform pgrad multiblock/proc
!
! Revision 1.1.1.1  2003/03/29 03:36:30  wasistho
! install ROCPERI
!
!
!******************************************************************************