PERI_coPgradUpdateFlu.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_coPgradUpdateFlu.F90,v 1.9 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 modbndpatch, ONLY   : t_patch
  USE moderror
  USE modmpi
  USE modparameters
  USE peri_modparameters
  IMPLICIT NONE

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

! ... loop variables
  INTEGER :: ijkc, ifc, ipatch

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

  INTEGER :: n, ijkc0, ijkcb, istage

  REAL(RFREAL)  :: rhoc, rhouc, rhovc, rhowc, rhoec, pc
  REAL(RFREAL)  :: refmassflux, rescaler, alpha, massflux
  REAL(RFREAL)  :: cfl, dt, rucm, delta, deltamassflux, aversurf
  REAL(RFREAL)  :: sndmassflux, rcvmassflux, sndaversurf,  rcvaversurf
  REAL(RFREAL)  :: snddt, rcvdt, avgdt, dtfactor, refmeanpgrad, reftauwall
  REAL(RFREAL)  :: tauwall, avgtauwall, sndtauwall, rcvtauwall
  REAL(RFREAL)  :: rn, sndrn, rcvrn, dy, fnmag, voltot
  REAL(RFREAL), POINTER :: xyz(:,:), vol(:), cofg(:,:)
  REAL(RFREAL), POINTER :: cv(:,:), tv(:,:), rdt(:)

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

  rcsidentstring = '$RCSfile: PERI_coPgradUpdateFlu.F90,v $ $Revision: 1.9 $'

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

! get dimensions, parameters and pointers -------------------------------------
 
  xyz  => region%grid%xyz
  vol  => region%grid%vol
  cofg => region%grid%cofg
  cv   => region%mixt%cv
  tv   => region%mixt%tv
  rdt  => region%dt

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

  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
    ENDIF
    n = 0
    DO ijkc0=1,region%grid%nCellsTot
      n = n+1
      avgdt = avgdt + rdt(ijkc0)
    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 -------------------------------------

  rucm   = 0._rfreal
  voltot = 0._rfreal
  DO ijkc = 1, region%grid%nCellsTot
    rucm   = rucm + cv(cv_mixt_xmom,ijkc0)*vol(ijkc)
    voltot = voltot + vol(ijkc)
  END DO
  massflux = rucm
  aversurf = 0.5_rfreal*voltot/delta

! this part only for channel flow

  avgtauwall = 0._rfreal

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

    DO ipatch=1,region%grid%nPatches
      patch => region%patches(ipatch)
      IF (patch%bcType>=bc_noslipwall .AND. &
          patch%bcType<=bc_noslipwall+bc_range) THEN   ! my boundary type

        DO ifc=1,patch%nBFaces
          ijkcb = patch%bf2c(ifc)

          dy    = abs( cofg(ycoord,ijkcb)-patch%fc(ycoord,ifc) )
          fnmag = patch%fn(xyzmag,ifc)

          tauwall = 0.5_rfreal*tv(tv_mixt_muel,ijkcb)* &
                    cv(cv_mixt_xmom,ijkcb)/cv(cv_mixt_dens,ijkcb)/dy
          avgtauwall = avgtauwall + tauwall*fnmag
        ENDDO ! ifc
      ENDIF   ! NOSLIPWALL
    ENDDO     ! iPatch  
  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
    sndaversurf = aversurf
    CALL mpi_allreduce( sndaversurf,rcvaversurf,1,mpi_double_precision,mpi_sum, &
                        global%mpiComm, global%mpierr )
    aversurf = rcvaversurf
  ENDIF
#endif
  IF (global%myProcId==masterproc) write(*,*)massflux,aversurf

! channel half width needed to update pressure gradient

  massflux   = massflux/aversurf
  avgtauwall = 0.5_rfreal*avgtauwall/aversurf

! update pressure gradient

  deltamassflux = refmassflux - massflux
  IF (global%myProcId==masterproc) WRITE(*,*) region%iRegionGlobal, &
    '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%iRegionGlobal, &
    ' 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

    rucm   = 0._rfreal
    DO ijkc = 1,region%grid%nCellsTot
      rucm   = rucm + cv(cv_mixt_xmom,ijkc0)*vol(ijkc)
    END DO
    massflux = rucm

#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/aversurf

    deltamassflux = refmassflux - massflux
    IF (global%myProcId==masterproc) WRITE(*,*) region%iRegionGlobal, &
      '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_coPgradUpdateFlu.F90,v $
! Revision 1.9  2008/12/06 08:44:37  mtcampbe
! Updated license.
!
! Revision 1.8  2008/11/19 22:17:49  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.7  2004/12/11 03:49:52  wasistho
! resereve dt=global%dtMin for dual tst
!
! Revision 1.6  2004/12/08 19:42:03  wasistho
! used dtFactor
!
! Revision 1.5  2004/12/08 01:47:50  wasistho
! bug fix, MPI_MAX to MPI_SUM
!
! Revision 1.4  2004/12/04 06:51:32  wasistho
! take global%dtMin for dt regardless steady or unsteady
!
! Revision 1.3  2004/06/17 20:04:02  wasistho
! compiled with RFLU
!
! Revision 1.2  2004/06/17 00:50:53  wasistho
! prepared for RFLU
!
! Revision 1.1  2004/06/09 01:10:26  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
!
!
!******************************************************************************