RFLU_ComputeExactFlowProbeError.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: Compute errors of computed solution relative to exact solution at 
!   probe locations.
!
! Description: None.
!
! Input: 
!   pRegion     Pointer to region
!
! Output: None.
!
! Notes: 
!   1. This routine assumes a perfect gas.
!
! ******************************************************************************
!
! $Id: RFLU_ComputeExactFlowProbeError.F90,v 1.3 2008/12/06 08:45:05 mtcampbe Exp $
!
! Copyright: (c) 2005 by the University of Illinois
!
! ******************************************************************************

SUBROUTINE rflu_computeexactflowprobeerror(pRegion)

  USE moddatatypes
  USE moderror
  USE moddatastruct, ONLY: t_region
  USE modmixture, ONLY: t_mixt_input
  USE modglobal, ONLY: t_global
  USE modgrid, ONLY: t_grid
  USE modparameters
    
  USE rflu_modbessel  
  USE rflu_modexactflow, ONLY:  rflu_computeexactflowpacoust 
  USE rflu_modflowhardcode, ONLY: rflu_getparamshardcodepacoust      
    
  USE modinterfaces, ONLY: mixtperf_d_cgp, & 
                           mixtperf_r_cpg
         
  IMPLICIT NONE

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

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

  TYPE(t_region), POINTER :: pregion

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

  LOGICAL :: printerrornorms
  CHARACTER(CHRLEN) :: ifilename,rcsidentstring
  INTEGER :: errorflag,ibc,icg,im,in,iprobe,iq
  REAL(RFREAL) :: atot,const,cpgas,dc,de,dinc,dtot,dummyreal,etaqm,ggas, &
                  idc,l,omega,pc,pe,probetime,ptot,rgas,ro,term,tc,ttot,uc, &
                  ue,vc,ve,wc,we,x,y,z
  REAL(RFREAL), DIMENSION(:,:), POINTER :: pcv,pdv
  TYPE(t_global), POINTER :: global
  TYPE(t_grid), POINTER :: pgrid
  TYPE(t_mixt_input), POINTER :: pmixtinput

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

  rcsidentstring = '$RCSfile: RFLU_ComputeExactFlowProbeError.F90,v $ $Revision: 1.3 $'

  global => pregion%global

  CALL registerfunction(global,'RFLU_ComputeExactFlowProbeError', &
                        'RFLU_ComputeExactFlowProbeError.F90')
 
  IF ( global%verbLevel > verbose_none ) THEN
    WRITE(stdout,'(A,1X,A)') solver_name, & 
      'Computing errors in flow solution at probe locations...'
                             
    IF ( global%verbLevel > verbose_low ) THEN 
      WRITE(stdout,'(A,3X,A,A)') solver_name,'Case: ',trim(global%casename)
    END IF ! global%verbLevel                         
  END IF ! global%verbLevel

! ==============================================================================
! Set pointers
! ==============================================================================

  pgrid      => pregion%grid
  pcv        => pregion%mixt%cv
  pdv        => pregion%mixt%dv
  pmixtinput => pregion%mixtInput

! ==============================================================================
! Set constants and initialize variables
! ==============================================================================

  cpgas = global%refCp
  ggas  = global%refGamma  
  rgas  = mixtperf_r_cpg(cpgas,ggas)
  
! ******************************************************************************
! Compute errors in probe quantities 
! ******************************************************************************

  SELECT CASE ( global%casename )

! ==============================================================================
!   Pipe acoustics. NOTE the pipe is assumed to have the x-coordinate 
!   running down the axis. 
! ==============================================================================

    CASE ( "pipeacoust" )
      CALL rflu_getparamshardcodepacoust(ptot,atot)
        dtot = mixtperf_d_cgp(atot,ggas,ptot)        

        l  = maxval(pgrid%xyz(xcoord,1:pgrid%nVert)) 
        ro = maxval(pgrid%xyz(ycoord,1:pgrid%nVert))        

        im  = max(pmixtinput%prepIntVal1,1)
        in  = max(pmixtinput%prepIntVal2,1)
        iq  = max(pmixtinput%prepIntVal3,1)
        ibc = max(min(pmixtinput%prepIntVal4,1),0)

        const = max(pmixtinput%prepRealVal1,0.0_rfreal)

        CALL rflu_jyzom(im,iq,dummyreal,etaqm,dummyreal,dummyreal)       

        omega = atot*sqrt((in*global%pi/l)**2 + (etaqm/ro)**2)         

        IF ( global%verbLevel > verbose_low ) THEN           
          WRITE(stdout,'(A,5X,A,1X,I2)'   ) solver_name, &
                'Boundary condition:',ibc          
          WRITE(stdout,'(A,5X,A,3(1X,I2))') solver_name, &
                'Mode:',im,in,iq                   
          WRITE(stdout,'(A,5X,A,1X,E13.6)') solver_name, &
                'Total density (kg/m^3):   ',dtot
          WRITE(stdout,'(A,5X,A,1X,E13.6)') solver_name, &
                'Total pressure (N/m^2):   ',ptot            
          WRITE(stdout,'(A,5X,A,1X,E13.6)') solver_name, &
                'Angular frequency (rad/s):',omega
          WRITE(stdout,'(A,5X,A,1X,E13.6)') solver_name, &
                'Constant (-):             ',const                  
        END IF ! global%verbLevel
  
! ------------------------------------------------------------------------------
!       Loop over probes and compute error: NOTE need to rewind because opening 
!       of probe file positions file at end.
! ------------------------------------------------------------------------------
  
        DO iprobe = 1,global%nProbes
          IF ( global%probePos(iprobe,probe_region) == & 
               pregion%iRegionGlobal ) THEN 
            icg = global%probePos(iprobe,probe_cell)

            rewind(if_probe+iprobe-1)

            WRITE(ifilename,'(A,I4.4)') trim(global%outDir)// & 
                                        trim(global%casename)//'.prbe_',iprobe
            OPEN(if_extr_data1,file=trim(ifilename),form='FORMATTED', &
                 status='UNKNOWN',iostat=errorflag)                            
            global%error = errorflag
            IF (global%error /= err_none ) THEN
              CALL errorstop(global,err_file_open,__line__, &
                             'File: '//trim(ifilename))
            END IF ! global%error

            emptyloop: DO 
              READ(if_probe+iprobe-1,*,iostat=errorflag) probetime, & 
                                                         dc,uc,vc,wc,pc,tc
              IF ( errorflag /= err_none ) THEN 
                EXIT emptyloop
              END IF ! errorFlag

              x = pgrid%cofg(xcoord,icg)
              y = pgrid%cofg(ycoord,icg)
              z = pgrid%cofg(zcoord,icg)

              CALL rflu_computeexactflowpacoust(global,z,y,x,probetime,l,ro, & 
                                                ibc,im,in,iq,etaqm,omega, &
                                                dtot,ptot,atot,const,de,ue, &
                                                ve,we,pe)  
                                                
              WRITE(if_extr_data1,'(1PE14.7,3(1X,E13.6))') & 
                    probetime,pe,pc,((pc-ptot)/(pe-ptot)-1.0_rfreal)
            END DO emptyloop
            
            CLOSE(if_extr_data1,iostat=errorflag)                            
            global%error = errorflag
            IF (global%error /= err_none ) THEN
              CALL errorstop(global,err_file_close,__line__, &
                             'File: '//trim(ifilename))
            END IF ! global%error                        
          END IF ! global%probePos  
        END DO ! iProbe

! ==============================================================================
!   Default - due to input error or missing CALL in this routine
! ==============================================================================  

    CASE default
      global%warnCounter = global%warnCounter + 1
    
      IF ( global%verbLevel > verbose_none ) THEN 
        WRITE(stdout,'(A,3X,A,1X,A)') solver_name,'*** WARNING ***', & 
              'Exact solution not available. Returning to calling procedure.'                                   
      END IF ! global%verbLevel
  END SELECT ! global%casename

! ******************************************************************************
! End
! ******************************************************************************

  IF ( global%verbLevel > verbose_none ) THEN 
    WRITE(stdout,'(A,1X,A)') solver_name, &
          'Computing errors in flow solution at probe locations done.'
  END IF ! global%verbLevel

  CALL deregisterfunction(global)

END SUBROUTINE rflu_computeexactflowprobeerror

! ******************************************************************************
!
! RCS Revision history:
!
! $Log: RFLU_ComputeExactFlowProbeError.F90,v $
! Revision 1.3  2008/12/06 08:45:05  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:18:16  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2005/04/29 12:41:32  haselbac
! Initial revision
!
! ******************************************************************************