module evalTijk
! Last updated 9/2/2007
!
! This module includes 5 subroutines used to evaluate a
! model of the three-point third order velocity correlation
! tensor for high Reynolds number isotropic turbulence
! as described in the Physics of Fluids paper entitled
!"An inertial range model for the three-point third-order velocity correlation"
! by Chang & Moser (2007). The model is expressed as the
! sum of a large number of terms which are rational in the magnitudes
! of the separation vectors. The coefficients and rational powers
! of these terms are read in from a data file, which is provided with
! these routines (basisdata.txt). See the document accompanying this
! module (Tijk_define.pdf) for details of the data representation.
!
! This set of subroutines is intended to provide an example of the evaluation
! of the 3-point 3rd-order correlation model. There is certainly room for
! improvement in these routines. Two obvious improvements are:
! 1) When evaluating multiple components of the tensor at the same r_ and s_
! it would be more efficient to evaluate them all at once.
! 2) It would be useful to implement an asymptotic form for the model that is
! valid when one of the separation vectors is small, so that under such
! conditions, evaluations are not polluted by roundoff error.
!
! The subroutines provided in this module (in the order in which
! they would normally be called) are:
!
! read_basis(datafile) : Reads the data describing the 5 basis tensors for
! the 3-point 3rd-order correlation. "datafile" is
! a string containing the name of the data file to read
!
! print_basis : Prints the basis characterization data to stdout
! in the same format as the input data file, though
! without comments. This is a debugging diagnostic.
!
! build_composite(weightsin,normalize)
! : Constructs a representation for a linear combination
! of the 5 basis tensors specified by weightsin
! an array of dimension 5 of double precision reals.
! This specifies a tensor in the 5-dimensional space
! of correlation tensors determined by the model.
! Normalize is a logical, if true, the weights will be
! normalized to sum to one, resulting in a correlation
! consistent with a dissipation of one. If not
! normalized, the dissipation related to the
! correlation is the sum of the coefficients in
! weightsin. If an individual basis is to be evaluated
! then weightsin should include all zeros, except for
! a one for the basis tensor to be evaluated.
!
! print_composite : prints the characterization of the composite,
! similar to print_basis. This is intended as a
! debugging diagnostic
!
!
! eval_composite(i,j,k,r_,s_,checkprecision)
! : is a function returning the value of
! T_{ijk}(r_,s_)=
! for the composite tensor defined by build_composite.
! i,j,k are integers, while r_ and s_ are
! three-dimensional double precision real vectors.
! checkprecision is a logical. If true, eval_composite
! will return a NaN if the magnitude the
! smallest separation vector
! is more than a factor of 100 smaller than
! the largest separation vector. This is useful because
! for such large aspect ratios, the evaluation
! of the model is subject to pollution by round-off
! errors. When this aspect ratio is less than 100,
! the roundoff error will be less than a part in 10^5.
implicit none
integer :: basis(64,3,14,5) !stores triple (coeff,rpow,spow) for
!up to 64 rational terms in r and s, for 14
!tensor components of 5 basis functions
integer :: nterms(14,5) !stores the number of rational terms for 14
!tensor components of 5 basis functions
integer :: basisnorms(5) !stores norms for 5 basis functions
double precision :: compositecoeffs(64,14) !stores coeff for up to 64
!rational terms in r and s, for 14 tensor
!components of the composite function
integer :: compositepowers(64,2,14) !stores (rpow,spow) for up to 64
!rational terms in r and s, for 14 tensor
!components of the composite function
integer :: compositeterms(14) !stores the number of rational terms for 14
!tensor components of the composite function
double precision :: w_sum
contains
subroutine read_comments
!internal subroutine used by read_basis
!ignores all lines that start with "!", carriage return lines,
!and lines with nothing but spaces
implicit none
character(500) :: str
do
str = ''
read(unit=7,fmt='(a)') str
if (scan(str,'!') == verify(str,' ') .or. len(trim(str)) == 0) then
! write (*,*) "ignoring: "//str
else
backspace(unit=7)
exit
end if
end do
end subroutine read_comments
subroutine read_basis(datafile)
!reads in basis functions from datafile, such as "basisdata.txt"
!comment lines (starting with "!") and carriage returns may be
!added to the datafile but order/format of the basis functions
!and order/format/spacing of tensor components must remain fixed
implicit none
integer :: i,j,k,ios,baseterms
character(*), intent(in) :: datafile
open(unit=7,file=datafile,iostat=ios,status="old",form="formatted",action="read")
if (ios/=0) then
print *,"unable to open file "
else
do k=1,5 !loop over basis functions
call read_comments
read (unit=7,fmt='(i8)') basisnorms(k)
baseterms = 0
do j=1,14 !loop over tensor components
call read_comments
read (unit=7,fmt='(i8)') nterms(j,k)
baseterms = baseterms + nterms(j,k)
call read_comments
do i=1,nterms(j,k) !loop over rational terms
read (unit=7,fmt='(i8)',advance="no") basis(i,1,j,k) !coeff
end do
read (unit=7,fmt=*)
call read_comments
do i=1,nterms(j,k) !loop over rational terms
read (unit=7,fmt='(i8)',advance="no") basis(i,2,j,k) !rpow
end do
read (unit=7,fmt=*)
call read_comments
do i=1,nterms(j,k) !loop over rational terms
read (unit=7,fmt='(i8)',advance="no") basis(i,3,j,k) !spow
end do
read (unit=7,fmt=*)
end do
write (*,*) "basis",k," total terms: ",baseterms
end do
close(unit=7)
end if
end subroutine read_basis
subroutine print_basis
!dumps to screen all basis data
!this should be the same as datafile but uncommented
!This is a debugging diagnostic
implicit none
integer :: i,j,k
do k=1,5
write (unit=*,fmt=*)
write (unit=*,fmt=*)
write (unit=*,fmt='(i8)') basisnorms(k)
do j=1,14
write (unit=*,fmt=*)
write (unit=*,fmt='(i8)') nterms(j,k)
do i=1,nterms(j,k)
write (unit=*,fmt='(i8)',advance="no") basis(i,1,j,k)
end do
write (unit=*,fmt=*)
do i=1,nterms(j,k)
write (unit=*,fmt='(i8)',advance="no") basis(i,2,j,k)
end do
write (unit=*,fmt=*)
do i=1,nterms(j,k)
write (unit=*,fmt='(i8)',advance="no") basis(i,3,j,k)
end do
write (unit=*,fmt=*)
end do
end do
end subroutine print_basis
subroutine build_composite(weightsin,normalize)
!This routine builds the description of the composite
!specified by the linear combination of the bases with
!coefficients in weightsin.
!input: weights for composite function, where
!compositefunction = SUM(weights(n)*basisfunction(n))
!output: compositecoeffs and compositepowers, coefficients
!and powers for rational functions in composite function
implicit none
double precision, intent(in) :: weightsin(5)
logical, intent(in) :: normalize
double precision :: weights(5)
double precision :: composite(14,-7:10,-7:10) !intermediate array of
!coefficents for 14 tensor components, 18 powers of r,
!and 18 powers of s
integer :: i,j,k,p,q
weights = weightsin
w_sum = weights(1)+weights(2)+weights(3)+weights(4)+weights(5)
!if necessary, normalize weights so that epsilon=1
if (w_sum /= 1. .and. w_sum /= 0. .and. normalize) then
weights = weights/w_sum
print *,"weights were normalized: ",weights
w_sum = 1
end if
composite = 0.
do k=1,5 !loop over basis functions
if (weights(k) /= 0.) then
do j=1,14 !loop over tensor components
do i=1,nterms(j,k) !loop over rational terms
!sum up coefficients in like terms (those with same tensor component, rpow, and spow) over 5 basis functions
composite(j,basis(i,2,j,k),basis(i,3,j,k)) = composite(j,basis(i,2,j,k),basis(i,3,j,k)) + basis(i,1,j,k)*weights(k)/basisnorms(k)
end do
end do
end if
end do
!build compositecoeffs and compositepowers from data in composite
do j=1,14 !loop over tensor components
i = 0
do p=-7,10 !loop over rpow
do q=-7,10 !loop over spow
if (composite(j,p,q) /= 0) then
i = i+1
if (i > 64) then
print *,"i > 64: compositecoeffs,compositepowers matrices too small"
else
compositecoeffs(i,j) = composite(j,p,q) !coeff
compositepowers(i,1,j) = p !rpow
compositepowers(i,2,j) = q !spow
end if
end if
end do
end do
compositeterms(j) = i
end do
end subroutine build_composite
subroutine print_composite
!dump composite description (uncommented)
!this if for diagnostics
implicit none
integer :: i,j
do j=1,14
write (unit=*,fmt=*)
write (unit=*,fmt='(i8)') compositeterms(j)
do i=1,compositeterms(j)
write (unit=*,fmt='(i8)',advance="no") compositecoeffs(i,j)
end do
write (unit=*,fmt=*)
do i=1,compositeterms(j)
write (unit=*,fmt='(i8)',advance="no") compositepowers(i,1,j)
end do
write (unit=*,fmt=*)
do i=1,compositeterms(j)
write (unit=*,fmt='(i8)',advance="no") compositepowers(i,2,j)
end do
write (unit=*,fmt=*)
end do
end subroutine print_composite
double precision function eval_composite(iin,jin,kin,rin_,sin_,checkprecision)
!input: i,j,k components for three velcities
!input: r_,s_ separation vectors
!output: numerical value of composite function at (r_,s_):
!T_{ijk}(r_,s_)=
implicit none
integer, intent(in) :: iin,jin,kin
double precision, intent(in) :: rin_(3),sin_(3)
logical, intent(in) :: checkprecision
integer :: i,j,k !components of three vectors
double precision :: r_(3),s_(3),t_(3) !separation vectors
double precision :: r,s,t !lengths of separation vectors
double precision :: rpow(-7:10), spow(-7:10) !r and s to 18 possible powers: -7,-6,...10
integer :: DELTA_ij,DELTA_jk,DELTA_ik !kronecker delta evaluated for input i,j,k
double precision :: sum1,sum2,minrst,maxrst,comp,tempreal
integer :: m,n,p,tempint
i=iin
j=jin
k=kin
r_=rin_
s_=sin_
!calculate lengths r, s, t, and vector t_
r = sqrt(r_(1)*r_(1) + r_(2)*r_(2) + r_(3)*r_(3))
s = sqrt(s_(1)*s_(1) + s_(2)*s_(2) + s_(3)*s_(3))
t_ = r_-s_
t = sqrt(t_(1)*t_(1) + t_(2)*t_(2) + t_(3)*t_(3))
if (i==j) then
DELTA_ij = 1
else
DELTA_ij = 0
end if
if (j==k) then
DELTA_jk = 1
else
DELTA_jk = 0
end if
if (i==k) then
DELTA_ik = 1
else
DELTA_ik = 0
end if
if (r==0.) then !evaluate 2pt correlations (much simpler)
eval_composite = w_sum/15*(s_(k)*DELTA_ij-1.5*(s_(i)*DELTA_jk+s_(j)*DELTA_ik))
else if (s==0.) then
eval_composite = w_sum/15*(r_(j)*DELTA_ik-1.5*(r_(k)*DELTA_ij+r_(i)*DELTA_jk))
else if (t==0.) then
eval_composite = -w_sum/15*(r_(i)*DELTA_jk-1.5*(r_(j)*DELTA_ik+r_(k)*DELTA_ij))
else
!evaluate full 3pt correlation
minrst = min(r,s,t)
maxrst = max(r,s,t)
if (minrst/maxrst < 0.01 .and. checkprecision) then
!print *,"minrst/maxrst < 0.01"
eval_composite = 0./0
else
if (minrst/=r .and. minrst/=s) then
!t is the smallest separation
!we found that evaluation as t->0 is less accurate than as r->0 or s->0
!therefore use symmetry T_{ijk}(r,s) = T_{jik}(-r,-t)
!swap subscripts and variables
tempint=i
i=j
j=tempint
r_ = -r_
s_ = -t_
!note we don't set t_ again since it is not used below
tempreal=s
s=t
t=tempreal
if (i==j) then
DELTA_ij = 1
else
DELTA_ij = 0
end if
if (j==k) then
DELTA_jk = 1
else
DELTA_jk = 0
end if
if (i==k) then
DELTA_ik = 1
else
DELTA_ik = 0
end if
end if
!normalize variables by t (T_{ijk} will be un-normalized at the end)
r_ = r_/t
s_ = s_/t
r = r/t
s = s/t
!pre-calculate r to different powers
rpow(0) = 1.
do p = 1,10
rpow(p) = rpow(p-1)*r
end do
do p = 1,7
rpow(-p) = 1./rpow(p)
end do
!pre-calculate s to different powers
spow(0) = 1.
do p = 1,10
spow(p) = spow(p-1)*s
end do
do p = 1,7
spow(-p) = 1./spow(p)
end do
sum2 = 0. !accumulates value of composite function
do n=1,14 !loop over tensor components
!determine which tensor component we are dealing with
select case (n)
case (1)
comp = r_(i)*r_(j)*r_(k)
case (2)
comp = r_(i)*r_(j)*s_(k)
case (3)
comp = r_(i)*s_(j)*r_(k)
case (4)
comp = s_(i)*r_(j)*r_(k)
case (5)
comp = s_(i)*s_(j)*s_(k)
case (6)
comp = s_(i)*s_(j)*r_(k)
case (7)
comp = s_(i)*r_(j)*s_(k)
case (8)
comp = r_(i)*s_(j)*s_(k)
case (9)
comp = r_(i)*DELTA_jk
case (10)
comp = r_(j)*DELTA_ik
case (11)
comp = r_(k)*DELTA_ij
case (12)
comp = s_(i)*DELTA_jk
case (13)
comp = s_(j)*DELTA_ik
case (14)
comp = s_(k)*DELTA_ij
end select
if (comp /= 0.) then
sum1 = 0. !accumulates value of rational function
do m=1,compositeterms(n) !sum over rational terms
sum1 = sum1 + compositecoeffs(m,n)*rpow(compositepowers(m,1,n))*spow(compositepowers(m,2,n))
end do
sum2 = sum2 + comp*sum1 !sum up tensor components times their respective rational functions
end if
end do
!un-normalized by t
eval_composite = sum2*t
end if
end if
end function eval_composite
end module evalTijk