TURB_InitSolution.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: Initialisation of TURB solution
!
! Description: Besides initialisation of TURB variables, in case non-uniform
!              filter is selected, non-uniform filter coefficients and
!              averaging coefficients are computed by calling the corresponding
!              routines. If wlm is active, define wlm mapping coefficients and
!              copy or interpolate wlm patch vals from level 1 to other levels.
!
! Input: region = data of current region
!
! Output: mueT, tCoT and cDyn get initial values. Wlm mapping coefficients
!         initiated and wlm patch vals interpolated to other levels.
!
! Notes: cDyn is set to zero every stage in ViscousFluxes.
!        Non-uniform quantities are recomputed in ViscousFluxes everytime
!        the grid moves.
!
!******************************************************************************
!
! $Id: TURB_InitSolution.F90,v 1.20 2009/08/26 12:28:52 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_initsolution( region ) ! PUBLIC

  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modgrid, ONLY       : t_grid
  USE modturbulence, ONLY : t_turb
  USE modglobal, ONLY     : t_global
  USE turb_modinterfaces, ONLY : turb_ranssageteddyvis, turb_wlminitia
#ifdef RFLO
  USE modinterfaces, ONLY : rflo_getdimensdummynodes, rflo_getnodeoffset
  USE turb_modinterfaces, ONLY : turb_flofacevolume, &
                                 turb_flofacewidth,  turb_flolesgencocc, &
                                 turb_flolesgencoff, turb_flowlmmetric 
#include "Indexing.h"
#endif
#ifdef RFLU
  USE modinterfaces, ONLY : mixtureproperties
  USE turb_modinterfaces, ONLY : turb_flufacevolume
#endif
  USE moderror
  USE modparameters
  USE turb_modparameters
  IMPLICIT NONE

! ... parameters
#ifdef RFLO
  TYPE(t_region) :: region
#endif
#ifdef RFLU
  TYPE(t_region), TARGET :: region
#endif

! ... loop variables
  INTEGER         :: ipatch, in, l, ifl

! ... local variables
  CHARACTER(CHRLEN)       :: rcsidentstring
  TYPE(t_global), POINTER :: global
  TYPE(t_grid) , POINTER  :: grid
  TYPE(t_turb), POINTER   :: turb
  TYPE(t_patch), POINTER  :: patch1, patch

  LOGICAL               :: dowlm
  INTEGER               :: turbmodel
  REAL(RFREAL)          :: triploc(xcoord:zcoord), treshold
  REAL(RFREAL), POINTER :: tv(:,:), tcv(:,:), tcvold(:,:), tdv(:,:), vort(:,:)
  REAL(RFREAL), POINTER :: dsterm(:,:)
#ifdef RFLO
  INTEGER :: ilev, inoff, ijnoff, ibn, ien
  INTEGER :: idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend, errorflag
#endif

  CHARACTER(2*CHRLEN+17) :: fname
  LOGICAL                :: fileexists

#ifdef RFLU
  INTEGER :: ibc, iec, ibn, ien, ic, ifg, ifgbeg
#endif
!******************************************************************************

  rcsidentstring = '$RCSfile: TURB_InitSolution.F90,v $'

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

! pre procedures -----------------------------------------------------------

#ifdef RFLU
  ibc = 1
  iec = region%grid%nCellsTot
  CALL mixtureproperties( region,ibc,iec,.true. )
#endif

! get dimensions -----------------------------------------------------------

#ifdef RFLO
  ilev =  region%currLevel
  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
                                 jdnbeg,jdnend,kdnbeg,kdnend )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
  ibn = indijk(idnbeg,jdnbeg,kdnbeg,inoff,ijnoff)
  ien = indijk(idnend,jdnend,kdnend,inoff,ijnoff)
#endif
#ifdef RFLU
  ibn = 1
  ien = region%grid%nFaces
#endif

! get parameters and pointers ----------------------------------------------

  turbmodel =  region%mixtInput%turbModel

#ifdef RFLO
  grid => region%levels(ilev)%grid
  turb => region%levels(ilev)%turb
  tv   => region%levels(ilev)%mixt%tv

  IF (region%turbInput%nDv > 0) THEN
    tdv => turb%dv
  ENDIF
  IF (region%turbInput%nCv > 0) THEN
    tcv    => turb%cv
    tcvold => turb%cvOld
    dsterm => turb%dsterm
  ENDIF
#endif
#ifdef RFLU
  grid => region%grid
  turb => region%turb
  tv   => region%mixt%tv

  IF (region%turbInput%nDv > 0) THEN
    tdv => region%turb%dv
  ENDIF
  IF (region%turbInput%nCv > 0) THEN
    tcv    => region%turb%cv
    tcvold => region%turb%cvOld
    dsterm => region%turb%dsterm
  ENDIF
#endif

! initialize turbulent viscosity and thermal conductivity at cells (tv of
! NS system), turbulence variables (turb.cvOld, dv, vorticities, etc),
! and metrics pertinent to turbulence

! general ---------------------------------------------------------
! first check if it's a restart from a laminar run
      IF (global%solutFormat == format_ascii .OR. global%solutFormat == format_hdf) THEN
        WRITE(fname,'(A,1PE11.5)') trim(global%inDir)//trim(global%casename)//'.turba_', &
                                   global%timeStamp
        INQUIRE(file=fname,exist=fileexists)
      ELSE IF (global%solutFormat == format_binary) THEN
        WRITE(fname,'(A,I6.6)') trim(global%inDir)//trim(global%casename)//'.turbb_', &
                                global%currentIter
        INQUIRE(file=fname,exist=fileexists)
      ELSE
        CALL errorstop( global,err_unknown_format,__line__ )
      ENDIF

  IF ((global%flowType == flow_unsteady .AND. &
       global%timeStamp <= 0._rfreal)    .OR. &
      (global%flowType == flow_steady   .AND. &
       global%currentIter <= 0)         .OR.  &
       (global%flowType == flow_unsteady .AND. &
       global%timeStamp > 0._rfreal   .AND. &
       fileexists .EQV. .false.)    .OR. &
       (global%flowType == flow_steady   .AND. &
       global%currentIter > 0     .AND. &
       fileexists .EQV. .false.)) THEN
    tv(tv_mixt_muet,:)  = 0._rfreal  
    global%esg1Sum      = 0._rfreal                ! pertinent to
    global%esg4Sum      = 0._rfreal                ! LES energy model 
    IF (region%turbInput%nCv > 0) THEN
#ifdef RFLO
      tcv(cv_sa_nutil,:)= tv(tv_mixt_muel,:)       ! pertinent to RaNS 
#endif
#ifdef RFLU
      tcv(cv_sa_nutil,:)= region%mixtInput%refVisc ! mixtProp, dv unknown yet
#endif
      tcvold(:,:)       = 0._rfreal          
    ENDIF
  ENDIF
  tv(tv_mixt_tcot,:)  = 0._rfreal

  IF (region%turbInput%nCv > 0) THEN
    dsterm(:,:) = 0._rfreal
  ENDIF

  IF (region%turbInput%nDv > 0) THEN
    tdv(:,:) = 0._rfreal
  ENDIF

  IF (ASSOCIATED( turb%vort )) THEN
    turb%vort = 0._rfreal
  ENDIF

! model dependent ---------------------------------------

  IF ((region%mixtInput%turbModel == turb_model_sa) .OR. &
      (region%mixtInput%turbModel == turb_model_dessa) .OR. &
      (region%mixtInput%turbModel == turb_model_hdessa)) THEN
    CALL turb_ranssageteddyvis( region ) ! nutilde to muet
  ENDIF

! metrics related ----------------------------------------

! compute face volumes needed for lesMij
  IF ((turbmodel==turb_model_fixsmag) .OR. &
      (turbmodel==turb_model_dynsmag) .OR. &
      (turbmodel==turb_model_dynmixd)) THEN
#ifdef RFLO
    CALL turb_flofacevolume( region,diri )
    CALL turb_flofacevolume( region,dirj )
    CALL turb_flofacevolume( region,dirk )
#endif
#ifdef RFLU
    CALL turb_flufacevolume( region )
#endif
  ENDIF

#ifdef RFLO
! compute filter coefficients

  IF (((turbmodel==turb_model_scalsim)  .OR. &
       (turbmodel==turb_model_dynsmag)  .OR. &
       (turbmodel==turb_model_dynmixd)) .AND. &
       (region%turbInput%filterType == filtype_nonunif)) THEN

    ALLOCATE( turb%workI(2,ibn:ien),stat=errorflag )
    ALLOCATE( turb%workJ(2,ibn:ien),stat=errorflag )
    ALLOCATE( turb%workK(2,ibn:ien),stat=errorflag )
    global%error = errorflag
    IF (global%error /= 0) CALL errorstop( global,err_allocate,__line__ )

    CALL turb_flofacewidth( region )
    CALL turb_flolesgencocc( region )
    CALL turb_flolesgencoff( region )
    DEALLOCATE( turb%workI,turb%workJ,turb%workK )
  ENDIF
#endif

! initialize zero-one fields

  IF (region%turbInput%nZof > 0) THEN
#ifdef RFLO
    turb%zofi = 1._rfreal
    turb%zofj = 1._rfreal
    turb%zofk = 1._rfreal
#endif
#ifdef RFLU
    turb%zofi  = 1._rfreal
    turb%bZofi = 1._rfreal
#endif
  ENDIF

! create zero-one fields

  IF (turbmodel==turb_model_fixsmag .OR. &
      turbmodel==turb_model_dynsmag .OR. &
      turbmodel==turb_model_dynmixd) THEN
    triploc(:) =  region%turbInput%xyzSmag(:)
    treshold   = -huge( 1.0_rfreal )/1000._rfreal

    IF ((triploc(xcoord) > treshold) .OR. &
        (triploc(ycoord) > treshold) .OR. &
        (triploc(zcoord) > treshold)) THEN

#ifdef RFLO
      DO in = ibn,ien
        DO l = xcoord,zcoord
          IF (grid%cfcI(l,in) < triploc(l)) THEN
            turb%zofi(l,zof_les_eddyvis,in) = 0._rfreal
          ELSE
            turb%zofi(l,zof_les_eddyvis,in) = 1._rfreal
          ENDIF
        ENDDO ! l
        DO l = xcoord,zcoord
          IF (grid%cfcJ(l,in) < triploc(l)) THEN
            turb%zofj(l,zof_les_eddyvis,in) = 0._rfreal
          ELSE
            turb%zofj(l,zof_les_eddyvis,in) = 1._rfreal
          ENDIF
        ENDDO ! l
        DO l = xcoord,zcoord
          IF (grid%cfcK(l,in) < triploc(l)) THEN
            turb%zofk(l,zof_les_eddyvis,in) = 0._rfreal
          ELSE
            turb%zofk(l,zof_les_eddyvis,in) = 1._rfreal
          ENDIF
        ENDDO ! l
      ENDDO   ! iN
#endif
#ifdef RFLU
      DO in = ibn,ien
        DO l = xcoord,zcoord
          IF (grid%fc(l,in) < triploc(l)) THEN
            turb%zofi(l,zof_les_eddyvis,in) = 0._rfreal
          ELSE
            turb%zofi(l,zof_les_eddyvis,in) = 1._rfreal
          ENDIF
        ENDDO ! l
      ENDDO   ! iN
      DO ipatch = 1,grid%nPatches
        patch  => region%patches(ipatch)
! TEMPORARY : removing usage of bf2bg from everywhere
!        ifgBeg =  patch%bf2bg(BF2BG_BEG)

        DO ifl = 1,patch%nBFaces
          ic  = patch%bf2c(ifl)
          ifg = ifl + ifgbeg-1

          DO l = xcoord,zcoord
            IF (grid%cofg(l,ic) < triploc(l)) THEN
              turb%bZofi(l,zof_les_eddyvis,ifg) = 0._rfreal
            ELSE
              turb%bZofi(l,zof_les_eddyvis,ifg) = 1._rfreal
            ENDIF
          ENDDO ! l
        ENDDO   ! ifl
      ENDDO     ! iPatch
#endif

    ENDIF ! triploc
  ENDIF   ! turbModel

! if applicable, initiate metric variables and utau of wlm

#ifdef RFLO
  DO ipatch=1,region%nPatches
    patch1 => region%levels(1)%patches(ipatch)
    patch  => region%levels(ilev)%patches(ipatch)
#endif
#ifdef RFLU
  DO ipatch=1,region%grid%nPatches
    patch  => region%patches(ipatch)
#endif

    dowlm = .false.
#ifdef RFLO
    IF (patch%bcType>=bc_noslipwall .AND. &
        patch%bcType<=bc_noslipwall+bc_range) THEN ! my boundary type
      IF (patch%valBola%switches(wlm_input_model) /= wlm_model_nomodel) THEN
        dowlm = .true.
      ENDIF
    ENDIF
#endif

    IF (dowlm .eqv. .true.) THEN

! --- get initial estimate of friction velocity utau
      CALL turb_wlminitia( region,patch )

#ifdef RFLO
! --- compute mapping coefficients from body fitted to cartesian and other metrics
      CALL turb_flowlmmetric( region,patch )

      IF (patch%valBola%distrib==bcdat_distrib) THEN
! ----- roughness distribution
!       CALL TURB_interpolate2Levels( )
        CALL errorstop( global,err_val_bcval,__line__,'No variable roughness yet'  )
      ELSE
! ----- distribution has constant value, so copied from a point on patch level 1
        patch%valBola%vals(:,wlm_vals_rough)= patch1%valBola%vals(1,wlm_vals_rough)
      ENDIF  ! distribution

      IF (patch%valBola%switches(wlm_input_model) == wlm_model_extern) THEN
! ----- wall stress distribution
!       CALL TURB_interpolate2Levels( )
      ENDIF
#endif

    ENDIF    ! doWlm
  ENDDO      ! iPatch

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

  CALL deregisterfunction( global )

END SUBROUTINE turb_initsolution

!******************************************************************************
!
! RCS Revision history:
!
! $Log: TURB_InitSolution.F90,v $
! Revision 1.20  2009/08/26 12:28:52  mtcampbe
! Ported to Hera.   Fixed logical expression syntax errors.  Replaced all
! IF (logical_variable)  with IF (logical_variable .eqv. .true.) as
! consistent with the specification.  Also changed: IF( ASSOCIATED(expr) )
! to IF ( ASSOCIATED(expr) .eqv. .true. ).   Intel compilers produce code
! which silently fails for some mal-formed expressions, so these changes
! are a net which should ensure that they are evaluated as intended.
!
! Revision 1.19  2009/08/12 04:15:59  mtcampbe
! Major update, bugfix from Abe development, more propagation compatibility,
! some Rocstar IO changes, Ju's temporary clipping fix for turbulence. A bug
! fix for initialization IO.
!
! Revision 1.18  2009/06/29 17:15:32  juzhang
! initialization of RANS/DES added for restart from a laminar run
!
! Revision 1.17  2008/12/06 08:44:41  mtcampbe
! Updated license.
!
! Revision 1.16  2008/11/19 22:17:53  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.15  2006/08/19 15:40:58  mparmar
! Removed bf2bg
!
! Revision 1.14  2006/01/30 23:06:15  wasistho
! removed comment after ifdef RFLU
!
! Revision 1.13  2006/01/17 17:51:40  wasistho
! applied tripping to all eddy viscosity models
!
! Revision 1.12  2006/01/12 23:46:55  wasistho
! POINTER to TARGET in rflu
!
! Revision 1.11  2006/01/12 09:48:49  wasistho
! enabled tripping fixed Smagorinsky
!
! Revision 1.10  2005/12/30 23:20:32  wasistho
! exclude rocflu from WLM treatment
!
! Revision 1.9  2005/12/29 19:48:01  wasistho
! modified rflu part
!
! Revision 1.8  2005/03/09 06:35:01  wasistho
! incorporated HDESSA
!
! Revision 1.7  2004/08/04 02:45:47  wasistho
! removed turb%avgCoI,J,K as it is defined as grid%c2fCoI,J,K
!
! Revision 1.6  2004/06/19 03:29:48  wasistho
! removed argument iReg in TURB_InitSolution
!
! Revision 1.5  2004/06/03 02:10:47  wasistho
! enabled non-uniform fix-Smagorinsky
!
! Revision 1.4  2004/03/20 03:28:29  wasistho
! prepared for RFLU
!
! Revision 1.3  2004/03/19 02:47:13  wasistho
! prepared for RFLU
!
! Revision 1.2  2004/03/12 02:55:35  wasistho
! changed rocturb routine names
!
! Revision 1.1  2004/03/05 04:37:00  wasistho
! changed nomenclature
!
! Revision 1.12  2004/02/26 21:26:29  wasistho
! initialize esg1Sum, esg4Sum
!
! Revision 1.11  2003/10/24 03:46:43  wasistho
! initiate mu_t to mu_l if RaNS is active
!
! Revision 1.10  2003/10/21 20:31:50  wasistho
! added dt relaxation in steady flow due to RANS source term
!
! Revision 1.9  2003/10/09 22:50:13  wasistho
! mv call to TURB_RansSAGetEddyVis from readSolution to initSolution
!
! Revision 1.8  2003/10/07 02:06:03  wasistho
! initial installation of RaNS-SA and DES
!
! Revision 1.7  2003/08/08 01:46:24  wasistho
! fixed turb. restart for steady flow
!
! Revision 1.6  2003/08/06 15:56:13  wasistho
! added vorticities computation
!
! Revision 1.5  2003/08/01 22:17:43  wasistho
! prepared rocturb for Genx
!
! Revision 1.4  2003/07/23 15:59:40  wasistho
! prepared more accurate rocturb restart
!
! Revision 1.3  2003/05/31 01:46:52  wasistho
! installed turb. wall layer model
!
! Revision 1.2  2002/10/16 07:48:19  wasistho
! Enable Fix Smagorinsky
!
! Revision 1.1  2002/10/14 23:55:29  wasistho
! Install Rocturb
!
!
!******************************************************************************