RFLU_TimeStepViscous.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: Calculate max. allowable local/global time step in the case
!          of viscous flow.
!
! Description: None.
!
! Input: 
!   pRegion             Pointer to region data
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: RFLU_TimeStepViscous.F90,v 1.16 2008/12/06 08:44:30 mtcampbe Exp $
!
! Copyright: (c) 2002-2004 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_timestepviscous(pRegion)

  USE moddatatypes
  USE moddatastruct, ONLY: t_region
  USE modglobal, ONLY: t_global
  USE modgrid, ONLY: t_grid
  USE modbndpatch, ONLY: t_patch
  USE moderror
  USE modparameters

  USE rflu_modgridspeedutils, ONLY: rflu_descalegridspeed
  USE rflu_modrindstates, ONLY: rflu_setrindstateinjectperf  
  
  USE modinterfaces, ONLY: mixtperf_c_grt, &
                           mixtperf_g_cpr, & 
                           mixtperf_r_m
  
  IMPLICIT NONE

! ******************************************************************************
! Declarations and definitions
! ******************************************************************************

! ==============================================================================
! Arguments
! ==============================================================================

  TYPE(t_region), POINTER :: pregion

! ==============================================================================
! Locals
! ==============================================================================

  CHARACTER(CHRLEN) :: rcsidentstring
  INTEGER :: c1,c2,distrib,gasmodel,icg,ifg,ifl,indcp,indgs,indmol,ipatch, & 
             turbmodel
  INTEGER :: dtminloc(1)
  INTEGER, DIMENSION(:), POINTER :: bf2c   
  INTEGER, DIMENSION(:,:), POINTER :: f2c 
  REAL(RFREAL) :: as,a1,a2,cp1,cp2,cp,dn,dtmin,dummyrealout,dx,dy,dz,fmue,fs, & 
                  fsu,gam,hl,ir1,ir2,minj,mm,mm1,mm2,nm,nx,ny,nz,pl,pr,prlam, &
                  prturb,rgas,rl,rr,r1,r2,srad,tinj,uinj,ul,u1,u2,vinj,vl,v1, & 
                  v2,winj,wl,w1,w2,vc                 
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,ptv,pvals
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

! ******************************************************************************
! Start
! ******************************************************************************

  rcsidentstring = '$RCSfile: RFLU_TimeStepViscous.F90,v $ $Revision: 1.16 $'

  global => pregion%global

  CALL registerfunction(global,'RFLU_TimeStepViscous',&
  'RFLU_TimeStepViscous.F90')

! ******************************************************************************
! Get dimensions and pointers
! ******************************************************************************

  pgrid => pregion%grid

  indgs = pgrid%indGs

  pcv => pregion%mixt%cv
  pdv => pregion%mixt%dv
  pgv => pregion%mixt%gv
  ptv => pregion%mixt%tv  
  
  prlam  = pregion%mixtInput%prLam
  prturb = pregion%mixtInput%prTurb

  indcp  = pregion%mixtInput%indCp
  indmol = pregion%mixtInput%indMol
  
  turbmodel = pregion%mixtInput%turbModel
  gasmodel  = pregion%mixtInput%gasModel

! ******************************************************************************
! Initialize time step
! ******************************************************************************

  DO icg = 1,pgrid%nCellsTot ! Explicit loop to avoid Frost problem
    pregion%dt(icg) = 0.0_rfreal 
  END DO ! icg

! ******************************************************************************
! Local time step
! ******************************************************************************

! ==============================================================================
! Interior faces
! ==============================================================================

  DO ifg = 1,pgrid%nFaces
    c1 = pgrid%f2c(1,ifg)
    c2 = pgrid%f2c(2,ifg)

    nx = pgrid%fn(xcoord,ifg)
    ny = pgrid%fn(ycoord,ifg)
    nz = pgrid%fn(zcoord,ifg)
    nm = pgrid%fn(xyzmag,ifg)                 

    dx = pgrid%cofg(xcoord,c2) - pgrid%cofg(xcoord,c1)
    dy = pgrid%cofg(ycoord,c2) - pgrid%cofg(ycoord,c1)
    dz = pgrid%cofg(zcoord,c2) - pgrid%cofg(zcoord,c1)            
    dn = sqrt(dx*dx + dy*dy + dz*dz)

    fs = pgrid%gs(indgs*ifg)

    r1 = pcv(cv_mixt_dens,c1)
    r2 = pcv(cv_mixt_dens,c2)        

    ir1 = 1.0_rfreal/r1
    ir2 = 1.0_rfreal/r2

    u1 = pcv(cv_mixt_xmom,c1)*ir1
    v1 = pcv(cv_mixt_ymom,c1)*ir1
    w1 = pcv(cv_mixt_zmom,c1)*ir1
    a1 = pdv(dv_mixt_soun,c1)

    u2 = pcv(cv_mixt_xmom,c2)*ir2
    v2 = pcv(cv_mixt_ymom,c2)*ir2
    w2 = pcv(cv_mixt_zmom,c2)*ir2
    a2 = pdv(dv_mixt_soun,c2)

    vc = 0.5_rfreal*((u1+u2)*nx + (v1+v2)*ny + (w1+w2)*nz)
    as = 0.5_rfreal*(a1+a2)

#ifdef TURB
    IF ( turbmodel == active ) THEN
      fmue = 0.5_rfreal*(ptv(tv_mixt_muel,c1) + ptv(tv_mixt_muel,c2))/prlam  & 
           + 0.5*rfreal*(ptv(tv_mixt_muet,c1) + ptv(tv_mixt_muet,c2))/prturb
    ELSE 
      fmue = 0.5_rfreal*(ptv(tv_mixt_muel,c1) + ptv(tv_mixt_muel,c2))/prlam
    END IF ! turbModel
#else
    fmue = 0.5_rfreal*(ptv(tv_mixt_muel,c1) + ptv(tv_mixt_muel,c2))/prlam
#endif      

    srad = (abs(vc - fs) + as)*nm & 
         + 4.0_rfreal*fmue/(0.5_rfreal*(r1 + r2))*nm/dn

    pregion%dt(c1) = pregion%dt(c1) + srad
    pregion%dt(c2) = pregion%dt(c2) + srad
  END DO ! ifg

! ==============================================================================
! Boundary faces
! ==============================================================================

  DO ipatch = 1,pgrid%nPatches
    ppatch => pregion%patches(ipatch)

    SELECT CASE ( ppatch%bcType ) 
    
! ------------------------------------------------------------------------------
!     Injection boundary
! ------------------------------------------------------------------------------   
    
      CASE ( bc_injection:bc_injection+bc_range ) 
        pvals => ppatch%mixt%vals     
         
        distrib = ppatch%mixt%distrib
        
        DO ifl = 1,ppatch%nBFaces
          c1 = ppatch%bf2c(ifl)

          nx = ppatch%fn(xcoord,ifl)
          ny = ppatch%fn(ycoord,ifl)
          nz = ppatch%fn(zcoord,ifl)
          nm = ppatch%fn(xyzmag,ifl)

          dx = pgrid%cofg(xcoord,c1) - ppatch%fc(xcoord,ifl)
          dy = pgrid%cofg(ycoord,c1) - ppatch%fc(ycoord,ifl)
          dz = pgrid%cofg(zcoord,c1) - ppatch%fc(zcoord,ifl)
          dn = sqrt(dx*dx + dy*dy + dz*dz)

          fs  = ppatch%gs(indgs*ifl)
          fsu = rflu_descalegridspeed(pregion,fs)

          rl = pcv(cv_mixt_dens,c1)
          pl = pdv(dv_mixt_pres,c1)

          minj = pvals(bcdat_inject_mfrate,distrib*ifl)

          IF ( minj > 0.0_rfreal ) THEN ! Surface burning
            tinj = pvals(bcdat_inject_temp,distrib*ifl)

            IF ( gasmodel == gas_model_tcperf ) THEN
              cp = pgv(gv_mixt_cp ,indcp *c1)        
              mm = pgv(gv_mixt_mol,indmol*c1)

              CALL rflu_setrindstateinjectperf(cp,mm,nx,ny,nz,minj,tinj,pl, &
                                               fsu,rl,ul,vl,wl,hl)
  
              r1   = rl
              vc   = ul*nx + vl*ny + wl*nz  
              rgas = mixtperf_r_m(mm)
              gam  = mixtperf_g_cpr(cp,rgas)              
              a1   = mixtperf_c_grt(gam,rgas,tinj)
            ELSE ! Defensive programming
              CALL errorstop(global,err_reached_default,__line__)              
            END IF ! gasModel
          ELSE ! Surface NOT burning
            r1 = rl
            vc = 0.0_rfreal
            a1 = pdv(dv_mixt_soun,c1)                        
          END IF ! minj
                  
#ifdef TURB
          IF ( turbmodel == active ) THEN
            fmue = ptv(tv_mixt_muel,c1)/prlam + ptv(tv_mixt_muet,c1)/prturb
          ELSE 
            fmue = ptv(tv_mixt_muel,c1)/prlam
          END IF ! turbModel
#else
          fmue = ptv(tv_mixt_muel,c1)/prlam
#endif  
                  
          pregion%dt(c1) = pregion%dt(c1) & 
                         + (abs(vc - fsu) + a1)*nm + 4.0_rfreal*fmue/r1*nm/dn
        END DO ! ifl

! ------------------------------------------------------------------------------
!     All other boundaries
! ------------------------------------------------------------------------------   
 
      CASE default        
        DO ifl = 1,ppatch%nBFaces 
          c1 = ppatch%bf2c(ifl)

          nx = ppatch%fn(xcoord,ifl)
          ny = ppatch%fn(ycoord,ifl)
          nz = ppatch%fn(zcoord,ifl)
          nm = ppatch%fn(xyzmag,ifl)

          dx = pgrid%cofg(xcoord,c1) - ppatch%fc(xcoord,ifl)
          dy = pgrid%cofg(ycoord,c1) - ppatch%fc(ycoord,ifl)
          dz = pgrid%cofg(zcoord,c1) - ppatch%fc(zcoord,ifl)
          dn = sqrt(dx*dx + dy*dy + dz*dz)

          fs  = ppatch%gs(indgs*ifl)
          fsu = rflu_descalegridspeed(pregion,fs)           

          r1 = pcv(cv_mixt_dens,c1)        

          ir1 = 1.0_rfreal/r1

          u1 = pcv(cv_mixt_xmom,c1)*ir1
          v1 = pcv(cv_mixt_ymom,c1)*ir1
          w1 = pcv(cv_mixt_zmom,c1)*ir1
          a1 = pdv(dv_mixt_soun,c1)

          vc = u1*nx + v1*ny + w1*nz

#ifdef TURB
          IF ( turbmodel == active ) THEN
            fmue = ptv(tv_mixt_muel,c1)/prlam + ptv(tv_mixt_muet,c1)/prturb
          ELSE 
            fmue = ptv(tv_mixt_muel,c1)/prlam
          END IF ! turbModel
#else
          fmue = ptv(tv_mixt_muel,c1)/prlam
#endif  

          pregion%dt(c1) = pregion%dt(c1) & 
                         + (abs(vc - fsu) + a1)*nm + 4.0_rfreal*fmue/r1*nm/dn
        END DO ! ifl                         
    END SELECT ! pPatch%bcType 
  END DO ! iPatch

! ******************************************************************************
! Finalize computation of local time step. NOTE set time step in dummy cells 
! to crazy value, that way code will blow up immediately should it be used by
! accident. 
! ******************************************************************************

  DO icg = 1,pgrid%nCells  
    pregion%dt(icg) = pgrid%vol(icg)/pregion%dt(icg)
  END DO ! icg

  DO icg = pgrid%nCells+1,pgrid%nCellsTot
    pregion%dt(icg) = REAL(ABS(CRAZY_VALUE_INT),kind=rfreal)
  END DO ! icg

! ******************************************************************************
! For unsteady flows, determine minimum step
! ******************************************************************************

  IF ( global%flowType == flow_unsteady ) THEN 
    dtmin    = minval(pregion%dt(1:pgrid%nCells))    
    dtminloc = minloc(pregion%dt(1:pgrid%nCells))
    
    pregion%dtMin    = dtmin
    pregion%dtMinLoc = dtminloc(1)    
  END IF ! global%flowType

! ******************************************************************************
! Finalize
! ******************************************************************************

  CALL deregisterfunction(global)

END SUBROUTINE rflu_timestepviscous

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_TimeStepViscous.F90,v $
! Revision 1.16  2008/12/06 08:44:30  mtcampbe
! Updated license.
!
! Revision 1.15  2008/11/19 22:17:43  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.14  2006/08/19 15:40:02  mparmar
! Renamed patch variables
!
! Revision 1.13  2006/04/07 15:19:22  haselbac
! Removed tabs
!
! Revision 1.12  2005/10/31 21:09:37  haselbac
! Changed specModel and SPEC_MODEL_NONE
!
! Revision 1.11  2005/05/28 18:03:56  haselbac
! Removed unnecessary setting of rgas and gamma
!
! Revision 1.10  2004/10/19 19:29:34  haselbac
! Now receive pRegion, partial rewrite, cosmetics
!
! Revision 1.9  2004/01/11 02:17:34  jiao
! Eliminated some redundant trailing spaces that made some lines too long.
! This changed was necessary to compile with NAG F90 compiler.
!
! Revision 1.8  2003/12/04 03:30:08  haselbac
! Changed timestep computation for viscous flows
!
! Revision 1.7  2003/10/15 02:43:57  haselbac
! Added code to determine location of minimum dt
!
! Revision 1.6  2003/08/28 00:58:43  haselbac
! Fixed bug in TURB section for timestep
!
! Revision 1.5  2003/04/18 20:00:09  haselbac
! Clean-up, added explicit initialization (prevent Frost-problem)
!
! Revision 1.4  2003/01/28 14:50:33  haselbac
! Use parameters in fn
!
! Revision 1.3  2002/11/27 20:27:45  haselbac
! Proper contribution from injecting boundaries
!
! Revision 1.2  2002/11/08 21:35:36  haselbac
! Fixed bug in use of grid speed
!
! Revision 1.1  2002/09/09 16:28:02  haselbac
! Initial revision
!
! ******************************************************************************