RFLU_TimeStepInviscid.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 inviscid flow.
!
! Description: None.
!
! Input: 
!   pRegion             Pointer to region
!
! Output: None.
!
! Notes: None.
!
! ******************************************************************************
!
! $Id: RFLU_TimeStepInviscid.F90,v 1.23 2008/12/06 08:44:30 mtcampbe Exp $
!
! Copyright: (c) 2002-2005 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_timestepinviscid(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

#ifdef PETSC
  USE rflu_modpetscnewtonkrylov, ONLY: rflu_petsc_getdtscale
#endif
  
  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
  INTEGER :: dtminloc(1) 
  REAL(RFREAL) :: as,a1,a2,cp,dtmin,dtscale,dummyrealout,fs,fsu,gam,hl,ir1, &
                  ir2,minj,mm,nm,nx,ny,nz,pl,pr,rgas,rl,rr,rur,rvr,rwr,srad, &
                  tinj,uinj,ul,u1,u2,vc,vinj,vl,vm2,v1,v2,winj,wl,w1,w2
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv,pgv,pvals
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_patch), POINTER :: ppatch

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

  rcsidentstring = '$RCSfile: RFLU_TimeStepInviscid.F90,v $ $Revision: 1.23 $'

  global => pregion%global

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

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

  pgrid => pregion%grid

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

  pcv => pregion%mixt%cv
  pdv => pregion%mixt%dv
  pgv => pregion%mixt%gv

  indgs = pgrid%indGs
  
! ******************************************************************************
! 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,pregion%grid%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)                  

    fs  = pgrid%gs(indgs*ifg)
    fsu = rflu_descalegridspeed(pregion,fs)    

    ir1 = 1.0_rfreal/pcv(cv_mixt_dens,c1)
    ir2 = 1.0_rfreal/pcv(cv_mixt_dens,c2)

    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)

    srad = (abs(vc - fsu) + as)*nm

    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)

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

          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)
  
              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
            vc = fsu
            a1 = pdv(dv_mixt_soun,c1)
          END IF ! minj
                  
          pregion%dt(c1) = pregion%dt(c1) + (abs(vc - fsu) + a1)*nm
        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)

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

          ir1 = 1.0_rfreal/pcv(cv_mixt_dens,c1)

          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

          pregion%dt(c1) = pregion%dt(c1) + (abs(vc - fsu) + a1)*nm
        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 Newton-Krylov solvers, ramp the timestep
! ******************************************************************************
 
  IF ( global%solverType == solv_implicit_nk ) THEN
    dtscale = 1.0
#ifdef PETSC
    CALL rflu_petsc_getdtscale(global,dtscale)
#endif
    DO icg = 1,pgrid%nCells
      pregion%dt(icg) = dtscale*pregion%dt(icg)
    ENDDO ! icg
  END IF ! global%solverType

! ******************************************************************************
! 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_timestepinviscid

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_TimeStepInviscid.F90,v $
! Revision 1.23  2008/12/06 08:44:30  mtcampbe
! Updated license.
!
! Revision 1.22  2008/11/19 22:17:43  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.21  2006/08/19 15:40:00  mparmar
! Renamed patch variables
!
! Revision 1.20  2006/04/07 15:19:22  haselbac
! Removed tabs
!
! Revision 1.19  2006/02/08 21:34:46  hdewey2
! Set dtScale so not undefined when running implicit solver wo PETSC
!
! Revision 1.18  2005/11/14 17:00:58  haselbac
! Removed DEBUG statement
!
! Revision 1.17  2005/10/31 21:09:37  haselbac
! Changed specModel and SPEC_MODEL_NONE
!
! Revision 1.16  2005/09/23 14:56:01  haselbac
! Bug fix: Missing ifdef PETSC
!
! Revision 1.15  2005/09/22 17:12:33  hdewey2
! Changed implicit timestep ramping to call RFLU_PETSC_GetDtScale to get the 
! timestep scaling.
!
! Revision 1.14  2005/08/17 20:28:24  hdewey2
! Added timestep ramping for the implicit Newton-Krylov solver
!
! Revision 1.13  2004/10/19 19:29:26  haselbac
! Now receive pRegion, partial rewrite, cosmetics
!
! Revision 1.12  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.11  2003/10/15 02:43:57  haselbac
! Added code to determine location of minimum dt
!
! Revision 1.10  2003/09/04 20:19:37  haselbac
! Removed temporary fix for Rocburn problem in GENx
!
! Revision 1.9  2003/08/22 20:31:56  haselbac
! Added temporary fix for ignition problems
!
! Revision 1.8  2003/04/18 20:00:09  haselbac
! Clean-up, added explicit initialization (prevent Frost-problem)
!
! Revision 1.7  2003/03/15 18:58:15  haselbac
! Clean-up, cosmetics
!
! Revision 1.6  2003/01/28 14:50:06  haselbac
! Use parameters in fn
!
! Revision 1.5  2002/11/27 20:27:45  haselbac
! Proper contribution from injecting boundaries
!
! Revision 1.4  2002/11/08 21:34:49  haselbac
! Fixed bug in use of grid speed
!
! Revision 1.3  2002/09/09 15:51:56  haselbac
! global and mixtInput under regions, add grid speeds
!
! Revision 1.2  2002/07/25 18:29:48  haselbac
! Bug fix, bug only shows for cases with dummy cells
!
! Revision 1.1  2002/05/04 17:02:00  haselbac
! Initial revision
!
! ******************************************************************************