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