RFLU_ComputeIntegral5OLES.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 integral 5 of optimal LES approach.
!
! Description: None.
!
! Input: 
!   pRegion     Pointer to current region
!
! Output: None.
!
! Notes: None.
!
!******************************************************************************
!
! $Id: RFLU_ComputeIntegral5OLES.F90,v 1.6 2008/12/06 08:44:29 mtcampbe Exp $
!
! Copyright: (c) 2002 by the University of Illinois
!
!******************************************************************************

SUBROUTINE rflu_computeintegral5oles(pRegion)

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

  USE modinterfaces, ONLY: rflu_computedcuhreinfo

  USE rflu_modoles

  IMPLICIT NONE

! --- parameters      

  TYPE(t_region), POINTER :: pregion

! --- local variables

  INTEGER :: c1g,c2g,c3g,c4g,errorflag,hloc,iint,ic1l,ic2l,ic3l,ic4l,ifc, & 
             ifcp,ifun,ifunnz,j,k,key,l,loopcounter,m,ncells,nfun,nfunnz, & 
             restartflag,vloc
  REAL(RFREAL) :: normfact,normfactterm
  
  TYPE(t_grid), POINTER :: pgrid    
  TYPE(t_global), POINTER :: global   
     
! ==============================================================================
! Start
! ==============================================================================

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

  IF ( global%myProcid == masterproc .AND. & 
       global%verbLevel > verbose_low ) THEN 
    WRITE(stdout,'(A,3X,A)') solver_name,'Computing integral 5...'     
  END IF ! global%verbLevel      

! ==============================================================================
! Set grid pointer
! ==============================================================================

  pgrid => pregion%grid

! ==============================================================================
! Set various quantities
! ==============================================================================

! ------------------------------------------------------------------------------
! Set DCUHRE parameters
! ------------------------------------------------------------------------------

  ndim   = 12
  nfun   = 81
  nfunnz = 21
  key    =  4

  errabsreq = 10.0_rfreal*epsilon(1.0_rfreal)
  errrelreq = 1.0e-2_rfreal 

! ------------------------------------------------------------------------------
! Compute DCUHRE information and allocate work array
! ------------------------------------------------------------------------------

  maxcalls = max_calls_limit

  CALL rflu_computedcuhreinfo(global,ndim,nfunnz,key,maxcalls,workarraysize)
  CALL rflu_allocatedcuhrearrays(global,ndim,nfunnz,nfun)
  
  ALLOCATE(integral(nfun),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'integral')
  END IF ! global%error

  ALLOCATE(workarray(workarraysize),stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_allocate,__line__,'workArray')
  END IF ! global%error

! ------------------------------------------------------------------------------
! Set normalization factor
! ------------------------------------------------------------------------------

  normfactterm = 1.0_rfreal/(pgrid%deltaOLES**12)

! ------------------------------------------------------------------------------
! Set non-zero integral components 
! ------------------------------------------------------------------------------

  CALL rflu_setmapfunnz2funcorr44(nfunnz)

! ==============================================================================
! Loop over protoype faces
! ==============================================================================

  DO ifcp = 1,3
    ifc = pgrid%fp2fOLES(ifcp)

! ------------------------------------------------------------------------------
!   Loop over cells
! ------------------------------------------------------------------------------        

    ncells = SIZE(pgrid%fsOLES,1)

    DO ic1l = 1,ncells
      c1g = pgrid%fsOLES(ic1l,ifc)

      lowlim(1) = pgrid%intLimOLES(int_lim_low,xcoord,c1g) 
      lowlim(2) = pgrid%intLimOLES(int_lim_low,ycoord,c1g)
      lowlim(3) = pgrid%intLimOLES(int_lim_low,zcoord,c1g)                

      upplim(1) = pgrid%intLimOLES(int_lim_upp,xcoord,c1g)
      upplim(2) = pgrid%intLimOLES(int_lim_upp,ycoord,c1g)
      upplim(3) = pgrid%intLimOLES(int_lim_upp,zcoord,c1g)              

      DO ic2l = 1,ncells
        c2g = pgrid%fsOLES(ic2l,ifc)

        lowlim(4) = pgrid%intLimOLES(int_lim_low,xcoord,c2g)
        lowlim(5) = pgrid%intLimOLES(int_lim_low,ycoord,c2g)
        lowlim(6) = pgrid%intLimOLES(int_lim_low,zcoord,c2g)              

        upplim(4) = pgrid%intLimOLES(int_lim_upp,xcoord,c2g)
        upplim(5) = pgrid%intLimOLES(int_lim_upp,ycoord,c2g)
        upplim(6) = pgrid%intLimOLES(int_lim_upp,zcoord,c2g)


        DO ic3l = 1,ncells
          c3g = pgrid%fsOLES(ic3l,ifc)

          lowlim(7) = pgrid%intLimOLES(int_lim_low,xcoord,c3g)
          lowlim(8) = pgrid%intLimOLES(int_lim_low,ycoord,c3g)
          lowlim(9) = pgrid%intLimOLES(int_lim_low,zcoord,c3g)            

          upplim(7) = pgrid%intLimOLES(int_lim_upp,xcoord,c3g)
          upplim(8) = pgrid%intLimOLES(int_lim_upp,ycoord,c3g)
          upplim(9) = pgrid%intLimOLES(int_lim_upp,zcoord,c3g)            

          DO ic4l = 1,ncells          
            c4g = pgrid%fsOLES(ic4l,ifc)

            lowlim(10) = pgrid%intLimOLES(int_lim_low,xcoord,c4g)
            lowlim(11) = pgrid%intLimOLES(int_lim_low,ycoord,c4g)
            lowlim(12) = pgrid%intLimOLES(int_lim_low,zcoord,c4g)                 

            upplim(10) = pgrid%intLimOLES(int_lim_upp,xcoord,c4g)
            upplim(11) = pgrid%intLimOLES(int_lim_upp,ycoord,c4g)
            upplim(12) = pgrid%intLimOLES(int_lim_upp,zcoord,c4g)       

! --------- Loop over integrals -----------------------------------------------

            DO iint = 1,3            
              errorflag   = 0
              maxcalls    = max_calls_start
              restartflag = 0 
              
              integral(:)   = 0.0_rfreal
              integralnz(:) = 0.0_rfreal              

              DO loopcounter = 1,dcuhre_loop_limit 
                IF ( iint == 1 ) THEN 
                  CALL dcuhre(ndim,nfunnz,lowlim,upplim,min_calls,maxcalls, & 
                              rflu_definecorrelation540,errabsreq,errrelreq, & 
                              key,workarraysize,restartflag,integralnz, & 
                              errabsest,neval,errorflag,workarray)            
                ELSE IF ( iint == 2 ) THEN 
                  CALL dcuhre(ndim,nfunnz,lowlim,upplim,min_calls,maxcalls, & 
                              rflu_definecorrelation541,errabsreq,errrelreq, & 
                              key,workarraysize,restartflag,integralnz, & 
                              errabsest,neval,errorflag,workarray)
                ELSE IF ( iint == 3 ) THEN 
                  CALL dcuhre(ndim,nfunnz,lowlim,upplim,min_calls,maxcalls, & 
                              rflu_definecorrelation542,errabsreq,errrelreq, & 
                              key,workarraysize,restartflag,integralnz, & 
                              errabsest,neval,errorflag,workarray)            
                END IF ! iInt

                IF ( errorflag == 0 .OR. loopcounter == dcuhre_loop_limit ) THEN
                  EXIT
                ELSE IF ( errorflag == 1 ) THEN 
                  restartflag = 1
                  maxcalls    = max_calls_factor*maxcalls

                  CALL rflu_computedcuhreinfo(global,ndim,nfunnz,key,maxcalls, & 
                                              workarraysizenew)

                  IF ( workarraysizenew > workarraysize ) THEN
                    EXIT       
                  END IF ! workArraySizeNew
                ELSE 
                  CALL errorstop(global,err_dcuhre_output,__line__)                            
                END IF ! errorFlag
              END DO ! loopCounter

! ----------- Normalize integral -----------------------------------------------

              IF ( iint == 1 ) THEN 
                normfact = normfactterm
              ELSE IF ( iint == 2 ) THEN 
                normfact = normfactterm*const_kolmogorov/ & 
                  (6.0_rfreal*pgrid%rhoOLES(ifc)**(2.0_rfreal/3.0_rfreal))
              ELSE IF ( iint == 3 ) THEN 
                normfact = normfactterm*const_kolmogorov*const_kolmogorov/ & 
                  (36.0_rfreal*pgrid%rhoOLES(ifc)**(4.0_rfreal/3.0_rfreal))         
              END IF ! iInt

              DO ifunnz = 1,nfunnz 
                ifun = mapfunnz2funcorr44(ifunnz)
                integral(ifun) = normfact*integralnz(ifunnz)
              END DO ! iFunNZ

! ----------- Store integral in array ------------------------------------------

              DO ifun = 1,nfun              
                CALL rflu_mapm2ijkl(ifun,l,m,j,k)

                vloc = rflu_getqposoles(j,k,ic3l,ic4l,ncells)
                hloc = rflu_geti4posoles(l,m,ic1l,ic2l,ncells)

                IF ( iint == 1 ) THEN
                  pgrid%int50OLES(ifcp,vloc,hloc) = integral(ifun) 
                ELSE IF ( iint == 2 ) THEN 
                  pgrid%int51OLES(ifcp,vloc,hloc) = integral(ifun)             
                ELSE 
                  pgrid%int52OLES(ifcp,vloc,hloc) = integral(ifun)             
                END IF ! iInt

              END DO ! iFun         

            END DO ! iInt

          END DO ! ic4l
        END DO ! ic3l
      END DO ! ic2l
    END DO ! ic1l       
  END DO ! ifcp


#ifdef CHECK_DATASTRUCT
! Data structure output for checking
  WRITE(stdout,'(A)') solver_name
  WRITE(stdout,'(A,1X,A)') solver_name,'### START CHECK OUTPUT ###'
  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES I50 integral matrix'
  DO ifcp = 1,3 ! loop over prototype faces
    WRITE(stdout,'(A,3X,A,1X,I2)') solver_name,'Face:',ifcp
    DO vloc = 1,9*ncells*ncells ! loop over components
      WRITE(stdout,'(A,1X,I6,36(1X,E11.4))') solver_name,vloc, & 
                                pgrid%int50OLES(ifcp,vloc,1:9*ncells*ncells)
    END DO ! vLoc
  END DO ! ifcp     
  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES I51 integral matrix'
  DO ifcp = 1,3 ! loop over prototype faces
    WRITE(stdout,'(A,3X,A,1X,I2)') solver_name,'Face:',ifcp
    DO vloc = 1,9*ncells*ncells ! loop over components
      WRITE(stdout,'(A,1X,I6,36(1X,E11.4))') solver_name,vloc, & 
                                pgrid%int51OLES(ifcp,vloc,1:9*ncells*ncells)
    END DO ! vLoc
  END DO ! ifcp      
  WRITE(stdout,'(A,1X,A)') solver_name,'Optimal LES I52 integral matrix'
  DO ifcp = 1,3 ! loop over prototype faces
    WRITE(stdout,'(A,3X,A,1X,I2)') solver_name,'Face:',ifcp
    DO vloc = 1,9*ncells*ncells ! loop over components
      WRITE(stdout,'(A,1X,I6,36(1X,E11.4))') solver_name,vloc, & 
                                pgrid%int52OLES(ifcp,vloc,1:9*ncells*ncells)
    END DO ! vLoc
  END DO ! ifcp
  WRITE(stdout,'(A,1X,A)') solver_name,'### END CHECK OUTPUT ###' 
  WRITE(stdout,'(A)') solver_name    
#endif 


! ==============================================================================
! Deallocate arrays
! ==============================================================================

  CALL rflu_deallocatedcuhrearrays(global)

  DEALLOCATE(integral,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'integral')
  END IF ! global%error

  DEALLOCATE(workarray,stat=errorflag)
  global%error = errorflag
  IF ( global%error /= err_none ) THEN 
    CALL errorstop(global,err_deallocate,__line__,'workArray')
  END IF ! global%error

! ==============================================================================
! End
! ============================================================================== 

  CALL deregisterfunction(global)     

END SUBROUTINE rflu_computeintegral5oles