vtf-logo

src/equations/spiralsgs.f

c----------------------------------------------------------------------
c
c Compressible/Incompressible stretched vortex subgrid model (LES)
c
c----------------------------------------------------------------------
c
c Copyright (C) California Institute of Technology 2000-2005
c Authors: D.I. Pullin  (dale@galcit.caltech.edu)
c          D.J. Hill
c          C.   Pantano
c
c----------------------------------------------------------------------

c     -----------------------------------------------------
c     Computes the circular average (in a given plane)
c     of the the subgrid contribution to the velocity variance
c     -----------------------------------------------------
c       This file needs to be edited when the plane is changed
c       or to change from circular averaged to spherical averaged
c     =====================================================
      Subroutine IntegrateStructure(ex,ey,ez,rx,ry,rz,StructureAve)
c     =====================================================
      implicit none

      include 'spiralsgs.i'
      
c     input
      double precision  ex, ey, ez
      double precision  rx, ry, rz
c      
c     output
      double precision  StructureAve
      double precision  normalcos,ratio
c
      if ( dir .eq. 1 ) then
c     --- structure function defined in the yz plane
         normalcos = ex
         ratio = rx
      else if ( dir .eq. 2 ) then
c     --- structure function defined in the xz plane
         normalcos = ey
         ratio = ry
      else if ( dir .eq. 3 ) then
c     --- structure function defined in the xy plane
         normalcos = ez
         ratio = rz
      endif
 
      if ( dir .gt. 0 ) then
c     For the circular average call this:      
         call IntegratePlaneStructure(normalcos,ratio,StructureAve)
      else
c     For a spherical average:
         StructureAve = 1.90695d0
      endif
      
      return
      end 

c     =====================================================
      subroutine ParticularStructureFunction(fstruc,vx,nvars,itype,
     &     dx,dy,dz,ixlo,ixhi,iylo,iyhi,izlo,izhi)
c     =====================================================

      implicit none

      include 'spiralsgs.i'

c     input
      integer  nvars, itype
      integer  ixlo,ixhi,iylo,iyhi,izlo,izhi
      double precision  vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
c     output
      double precision  fstruc(ixlo:ixhi,iylo:iyhi,izlo:izhi)
      
      integer  dir1, dir2,dir3
      double precision  weight
      double precision  strucr
      double precision  dx,dy,dz
      double precision  ddir1,dir1r_23
      double precision  ddir2,dir2r_23
      double precision  ddir3,dir3r_23
      integer  i,j,k

c--------------------------------------------------
      if ( dir .eq. 0 ) then
c     --- structure function defined in a sphere
         dir1 = 1
         ddir1 = dx
         dir2 = 2
         ddir2 = dy
         dir3 = 3 
         ddir3 = dz
      else if ( dir .eq. 1 ) then
c     --- structure function defined in the yz plane
         dir1 = 2
         ddir1 = dy
         dir2 = 3
         ddir2 = dz
c     don't use third direction      
         dir3 = 0   
         ddir3 = dx      
      else if ( dir .eq. 2 ) then
c     --- structure function defined in the xz plane
         dir1 = 1
         ddir1 = dx
         dir2 = 3
         ddir2 = dz
c     don't use third direction
         dir3 = 0   
         ddir3 = dy
      else if ( dir .eq. 3 ) then
c     --- structure function defined in the xy plane
         dir1 = 2
         ddir1 = dy
         dir2 = 1
         ddir2 = dx
c     don't use third direction
         dir3 = 0   
         ddir3 = dz
      endif
c__________________________________________________

      if ( dir .eq. 0 ) then
c      spherical average:
         strucr = (ddir1*ddir2*ddir3)**(1.D0/3.D0)
      else
c      circular average:      
         strucr = dsqrt(ddir1*ddir2)
      endif

      dir1r_23 = (strucr/ddir1)**(2.D0/3.D0)
      dir2r_23 = (strucr/ddir2)**(2.D0/3.D0)
      dir3r_23 = (strucr/ddir3)**(2.D0/3.D0)

      if (dir.eq.0) then
c      spherical average
         weight = 1.0d0/6.0d0
      else
c      circular average:  
         weight = 0.25d0
      end if
c     
      do k = izlo, izhi
         do j = iylo,iyhi
            do i = ixlo,ixhi
               fstruc(i,j,k) = 0.0d0
            enddo
         enddo
      enddo
      
      call compute_square_diffs(fstruc,vx,nvars,itype, 
     &     ixlo,ixhi,iylo,iyhi,izlo,izhi,
     &     dir1r_23, dir1)
c
      call compute_square_diffs(fstruc,vx,nvars,itype, 
     &     ixlo,ixhi,iylo,iyhi,izlo,izhi,
     &     dir2r_23, dir2)
c     
      if (dir.eq.0) then
c       -- we are doing the spherical average, need one more direction
         call compute_square_diffs(fstruc,vx,nvars,itype, 
     &        ixlo,ixhi,iylo,iyhi,izlo,izhi,
     &        dir3r_23, dir3)
      end if

      do k = izlo, izhi
         do j = iylo,iyhi
            do i =ixlo,ixhi
               fstruc(i,j,k) = weight * fstruc(i,j,k)
            end do
         end do
      end do

c      
      return
      end



c     =====================================================
      subroutine IntegratePlaneStructure(e3z,ratio,q1)
c     =====================================================
c ----- This computes the integral (B8) 
c ----- in Voelkl,Pullin,Chan Phy Fluids 2000
c ----- and mutiplies it by 2/pi - giving the circular
c ----- average of the "resolved scales"
      implicit none

c     input
      double precision  e3z
      double precision  ratio
c     output
      double precision  q1
c
      double precision pi
c
      double precision  calpha,ssqalpha,afac,cfac,sapow,ctmp
c
      pi=3.141592653589793238462643383279502884197

      calpha = ABS(e3z)
c     this is ssqalpha = sin\psi on Page 1823 of Voelkl,Pullin,Chan
c     \simga
      ssqalpha = 1.D0 - calpha*calpha

c     (3 \pi^{7/3} / 16) (2-\sigma)      
      afac = 2.7103 * (2.D0 - ssqalpha)

c     C_1 G(\sigma)
c     where G(\sigma) = \pi(2 -1/3 \sigma - 1/12 \sigma^2 
c                            -25/648\simga^3 -175/776 \sigma^4 )
c
c     and C_1 = \pi /(2^(2/3) \Gamma^2(4/3)
      cfac = 9.004 - 1.501*ssqalpha
      sapow = ssqalpha * ssqalpha
      cfac = cfac - 0.375*sapow
      sapow = sapow * ssqalpha
      cfac = cfac - 0.174*sapow
      sapow = sapow * ssqalpha
      cfac = cfac - 0.101*sapow
      
c
      ctmp = (cfac + afac * ratio**(1.33333333333) )
      if ( ctmp .ne. 0.0d0 ) then
         q1 = afac * cfac * ratio * ratio / ctmp
      else
         q1 = 0.0d0
      endif
c
      q1 = q1*2.0/pi

      return
      end

c     =====================================================      
      subroutine compute_square_diffs(temp, vx, nvars, itype, 
     &     ixlo,ixhi,iylo,iyhi,izlo,izhi,dr_23, id)
c     =====================================================

      implicit none
      
      integer  nvars, itype
      integer  ixlo,ixhi,iylo,iyhi,izlo,izhi
      integer  id
      double precision  dr_23
   
      double precision  vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
      double precision  temp(ixlo:ixhi,iylo:iyhi,izlo:izhi)

      integer i,j,k
      integer i1, i2, i3
      
      if ( itype .eq. 0 ) then
         i1 = 2
         i2 = 3
         i3 = 4
      endif

      if (id.eq.1) then

         if ( itype .eq. 0 ) then
            do k = izlo, izhi
               do j = iylo,iyhi
                  do i = ixlo+1,ixhi-1
                     temp(i,j,k) = temp(i,j,k) + (
     &                    (vx(i1,i+1,j,k)-vx(i1,i,j,k))**2 +  
     &                    (vx(i1,i-1,j,k)-vx(i1,i,j,k))**2 +  
     &                    (vx(i2,i+1,j,k)-vx(i2,i,j,k))**2 +  
     &                    (vx(i2,i-1,j,k)-vx(i2,i,j,k))**2 +
     &                    (vx(i3,i+1,j,k)-vx(i3,i,j,k))**2 +  
     &                    (vx(i3,i-1,j,k)-vx(i3,i,j,k))**2 ) * dr_23
                  end do
               end do
            end do
         else
            do k = izlo, izhi
               do j = iylo,iyhi
                  do i = ixlo+1,ixhi-1
                     temp(i,j,k) = temp(i,j,k) + (
     &                    (vx(itype,i+1,j,k)-vx(itype,i,j,k))**2 +  
     &                    (vx(itype,i-1,j,k)-vx(itype,i,j,k))**2) 
     $                    * dr_23 
                  end do
               end do
            end do
         endif
         
      else if (id.eq.2) then
         
         if ( itype .eq. 0 ) then
            do k = izlo, izhi
               do j = iylo+1,iyhi-1
                  do i = ixlo,ixhi
                     temp(i,j,k) = temp(i,j,k) + (
     &                    (vx(i1,i,j+1,k)-vx(i1,i,j,k))**2 +  
     &                    (vx(i1,i,j-1,k)-vx(i1,i,j,k))**2 +  
     &                    (vx(i2,i,j+1,k)-vx(i2,i,j,k))**2 +  
     &                    (vx(i2,i,j-1,k)-vx(i2,i,j,k))**2 +
     &                    (vx(i3,i,j+1,k)-vx(i3,i,j,k))**2 +  
     &                    (vx(i3,i,j-1,k)-vx(i3,i,j,k))**2 ) * dr_23
                  end do
               end do
            end do
         else
            do k = izlo, izhi
               do j = iylo+1,iyhi-1
                  do i = ixlo,ixhi
                     temp(i,j,k) = temp(i,j,k) + (
     &                    (vx(itype,i,j+1,k)-vx(itype,i,j,k))**2 +  
     &                    (vx(itype,i,j-1,k)-vx(itype,i,j,k))**2 ) 
     $                    * dr_23
                  end do
               end do
            end do
         endif
         
      else 

         if ( itype .eq. 0 ) then
            do k = izlo+1, izhi-1
               do j = iylo,iyhi
                  do i = ixlo,ixhi
                     temp(i,j,k) = temp(i,j,k) + (
     &                    (vx(i1,i,j,k+1)-vx(i1,i,j,k))**2 +  
     &                    (vx(i1,i,j,k-1)-vx(i1,i,j,k))**2 +  
     &                    (vx(i2,i,j,k+1)-vx(i2,i,j,k))**2 +  
     &                    (vx(i2,i,j,k-1)-vx(i2,i,j,k))**2 +
     &                    (vx(i3,i,j,k+1)-vx(i3,i,j,k))**2 +  
     &                    (vx(i3,i,j,k-1)-vx(i3,i,j,k))**2 ) * dr_23
                  end do
               end do
            end do
         else
            do k = izlo+1, izhi-1
               do j = iylo,iyhi
                  do i = ixlo,ixhi
                     temp(i,j,k) = temp(i,j,k) + (
     &                    (vx(itype,i,j,k+1)-vx(itype,i,j,k))**2 +  
     &                    (vx(itype,i,j,k-1)-vx(itype,i,j,k))**2 ) 
     $                    * dr_23
                  end do
               end do
            end do
         endif
         
      end if
      
      return
      end 
      

<