\documentstyle[11pt,epsf]{report}
\setlength{\textwidth}{6.6truein}
\setlength{\oddsidemargin}{-.2truein}
\setlength{\evensidemargin}{-.2truein}
\setlength{\textheight}{9truein}
\setlength{\topmargin}{-.4truein}
\setlength{\headsep}{.2truein}
\setlength{\footskip}{.3truein}
\begin{document}

\title{
Towards a Fractal Representation of Turbulent Flows}
\author{Leon van Dommelen}
\date{FMRL Report LVD98A, 3/5/98}
\maketitle

\begin{abstract}

Fractal representation of turbulent flows is considered.  The
objective is to use fractal techniques to approximate a larger range
of scales in a more accurate fashion than would be possible by simple
models of the small scales. We envision that fractal representation
will eventually allow complex high Reynolds number turbulent flows to
be approximated at all scales without statistically significant errors
while using only a small set of numerical parameters.  Such a reduced
set of parameters would facilitate storage, retrieval, and management
of these flows.  More importantly, the fractal representation could be
mined for information about the properties of the turbulent flow. This
information could be implemented in large eddy simulations through the
subgrid model.  In contrast to turbulence models that attempt to
formulate universal descriptions for all turbulent flows, the proposed
techniques are based on a search for local similarities in specific
flows. Much remains to be done to achieve those objectives.

\end{abstract}

\chapter{Introduction}
\label{sec=intro}

Without doubt, the most important practical problem in turbulence
is that the number of variables needed to describe them fully tends to
be prohibitive.  Accepting that the unsteady three-dimensional
Navier-Stokes equations are properly posed, they should be numerically
solvable given enough computational power.  However, such power tends
to be out of reach, and it will probably remain a very wasteful
approach even if it should eventually become possible.

The most common approach to deal with this problem is to use some kind
of averaging to eliminate the small scale flow processes.  The
resulting smoothed flow field can be described using much less
computational degrees of freedom.  However, due to the nonlinearity of
the Navier-Stokes equations, the small scale processes are not truly
eliminated but do still affect the smoothed flow field.
To resolve this, crude approximations are made for the effects of the
unknown small-scale processes.  The result is a replacement of the
true flow by a much more easily described smooth solution of an
artificial equation \cite{lesieur}.

In a large eddy simulation (LES), the hope is that if the computation
resolves enough of the large-scale flow processes, the behavior of the
small scale ones may be sufficiently general that it can be described
by some simple model.  However, the evidence suggests that the subgrid
effects are not universal, but depend significantly on the type of
flow field.  Hence, they may not be easy to model well for completely
general flows.  It has also been found that current subgrid scale
models do not provide an accurate description of the effects of the
subgrid scales on the larger scales, but only simulate their net
effects \cite{piomelli}.

Even more worrisome, the results of LES seems to deteriorate
significantly when the range of modeled scales increases.  Hence the
development of better subgrid scale models is seen as a pacing item
for the practical application of the very large eddy simulations
(VLES) which are of the greatest interest \cite{ferziger}.

The experience with Reynolds averaging approaches seems to suggest
that while various models of the small scales often work
satisfactorily, it may not be easy to judge their reliability or to
significantly improve upon them when higher accuracy is called for.
Although a large eddy simulation may, hopefully, make the modeling
easier and less critical, it is not clear what to do in case the model
would not be sufficient for the required accuracy.

An alternative approach to be explored here is to keep all scales, but
to reduce the number of degrees of freedom by an alternate assumption.
The fact that a very large range of scales is of importance does not
automatically mean that the required number of parameters to describe
them sufficiently accurately increases proportionally.  There may be
much less parameters in a suitable unconventional description.

Examples are mathematical objects known as self-similar fractals.
These fractals contain an infinity of scales, and a straightforward
mesh-based description of them would require infinite resolution
\cite{peitgen}.  A simple classical fractal is the Cantor set shown in
figure~\ref{fig=cantor0}; it consists of all numbers between zero and
one that can be written without using the digit ``1'' in a base 3
representation.  For example, the middle third of the interval is not
included since in base 3, the numbers start there as ``0.1...''.  In
general, the size of the holes in the Cantor set is determined by the
location of the first digit 1; the further it is behind the decimal
point, the smaller the hole. Thus the Cantor set has infinitely fine
detail.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[40 360 652 450]{cantor0.ps}}
  \caption{The Cantor set.}
  \label{fig=cantor0}
\end{figure}

Because of this detail, a ``mesh-based'' representation of the Cantor
set using, say, piecewise polynomials or a Fourier series requires
infinitely many polynomials, terms, etcetera.  There is however a very
powerful, quite general, way of representing this type of fractals:
instead of describing the object itself, one can describe a machine
that makes it.

The machine that makes the Cantor set works as follows: it accepts any
arbitrary set of points and replaces it by two reduced copies of that
set.  The reduction factor is three and the reduced copies are spaced
a distance 2/3 apart.  Figure~\ref{fig=cantor} show the repeated
application of this machine to the unit interval; the results converge
to the Cantor set.  A machine of this type is called an IFS (for
Iterated Function System).

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[190 200 502 400]{cantor.ps}}
  \caption{IFS convergence to the Cantor set.}
  \label{fig=cantor}
\end{figure}

All the properties of the Cantor set are in the simple machine.  The
initial condition is unimportant; the process would converge to the
Cantor set for any other initial set instead of the unit interval.  A
surprisingly large number of other shapes can be created by merely
changing the parameters of the IFS, (the number of reduced copies,
their reduction factors and placement, and reduction method), A
classical example is the fern shape created by
Barnsley \cite{barnsley}.

By construction, the Cantor set consists of parts that are similar to
the whole set; the small scales are simply reduced copies of the large
scales.  This is different for turbulent flows, where the Reynolds
number for the smaller scale processes is smaller.  If a turbulent
flow is to be approximated by an IFS, it is necessary to incorporate
Reynolds number effects by some means of smoothing the small
scales. Smoothing of the small scales should result in an IFS with
scales limited by Kolmogorov limits.  One important issue is how to
perform such smoothing.

Another issue is that even without Reynolds number effects, turbulent
flows are not self-similar. Further, in the presence of boundaries,
arbitrary regions, say at positions for which the flow geometry is not
similar, cannot be expected to be similar even in some statistical
sense.  The solution adopted here is to use partitioning of the flow.
This allows the construction of partitioned IFS (PIFS) representations
in which all desired large scales of turbulent flows can be
represented arbitrarily accurately.

The IFS will fill in the unrepresented smaller scales with reduced
smoothed copies of larger scales taken from regions that resemble the
local flow.  Of course, after severely reducing the true number of
degrees of freedom of a turbulent flow, all possible realizations of
the small scales can no longer be represented.  However, the hope is
that with suitable corrections for Reynolds number effects, the
partitioned representation will give one possible realization with
statistics close to the ensemble average statistics.

A line of reasoning that favors such optimism is that if two regions
are very similar in all their major statistics, it seems reasonable to
assume that they are also very similar in their more rapidly varying
smaller scales.  Note the more restricted scope compared to an
ordinary LES: in this case we only search for similarities between
regions that already appear almost the same in their global
properties.  This reduced need for generality suggests that the
fractal representation may eventually turn out to be more reliable,
accurate, and robust than a direct modeling approach.

Computing flows completely using only PIFS seems a very challenging
proposition.  Our less ambitious first goal is to attempt to use PIFS
representation as part of more conventional LES schemes.  From the
computed large scales, a PIFS representation can be formed. Sampling
of the small scales predicted by this PIFS may then be used as a guide
in adjusting the subgrid scale model in the LES.

To give an heuristic justification: suppose that a LES contains a
region that is sufficiently large so that there is still confidence in
its somewhat smaller scales, and it is seen that these smaller scales
are, say, noticeably less than would be expected from the used subgrid
model.  In that case, would it not be reasonable to suppose that an
other, somewhat smaller, region which resembles it closely would also
have subgrid scales that would be smaller than predicted?  This idea
resembles the dynamic subgrid model of Germano \cite{lesieur} somewhat,
but it is based only on observation of flow features, not on the
postulated validity of a particular subgrid scale model.

\chapter{Iterated Function Systems}
\label{sec=ifs}

In this section we will briefly discuss some basic properties of
iterated function systems such as the one that generated the Cantor
set in figure~\ref{fig=cantor}.  We also discuss partitioning as
a way to avoid the limitations of self similarity.

\section{The contraction mapping principle}
\label{sec=cm}

The contraction principle is the theoretical justification for IFS of
the type used in this study.  To understand the principle and how it
applies, assume that an IFS machine is used on two different sets.
The difference between the two sets becomes smaller by applying the
machine on them. For example, the machine that creates the Cantor set,
figure~\ref{fig=cantor}, will replace each of the sets by two copies
that are reduced by a factor 3: the distances between the sets are
then at least reduced by the same factor 3. With successive
application, the distance between the sets converges geometrically.

It must follow that all initial sets converge to the same final
``attractor'' set.  In the example, the final set is always the Cantor
set, regardless of the initial (nonempty compact) set.  In general,
the attractor is unique; it is the single set that is unchanged by the IFS
machine.

(For a more complete description, including the definition of
the Hausdorff distance between sets that gives rise to a closed
compact space of sets, see \cite{peitgen}.  Put simply, the Hausdorff
distance is the largest distance of any point in either set to the
other set.)

\section{Smoothed iterated function systems}
\label{sec=sifs}

Using an IFS, it is possible to create infinitely fine detail.
However, in turbulent flows, viscous diffusion limits the smallest
scales that can exist.  To simulate this using an IFS, it is possible
to smooth the image obtained by the mapping.

For the IFS leading to the Cantor set, a simple smoothing can be
achieved by adding all points within some small distance $\epsilon$ to
the transformed set.  Since the length of each continuous segment of
the set is now at least 2$\epsilon$, the set will consists of a finite
number of segments.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[190 230 502 400]{cantors.ps}}
  \caption{SIFS convergence.}
  \label{fig=cantors}
\end{figure}

The contraction mapping principle is not affected by this smoothing
since it will not increase the distances between sets, in a linear
normed vector space.  For the smoothings used for the turbulent
flows in this study, this is no longer true.  However, experimentally
the constructed IFS always converged.

\section{Partitioned iterated function systems}
\label{sec=pifs}

Iterated functions systems such as the one that produces the Cantor
set require that the set is made up from reduced copies of itself.
Another simple example of a self-similar ``fractal'' is the line in
figure \ref{fig=liness}. It can be thought of as consisting of two
equal halves, where each half is a copy of the full line reduced by a
factor two. Hence, the IFS machine that produces this line generates
two half-size copies, placing one at the origin and another at
(0.5,1). If this IFS is modified by additionally reflecting the second
half horizontally, the fractal figure~\ref{fig=fracss} results.  Again
this fractal is self-similar: the halves are still reduced copies of
the full curve.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{liness.ps}}
  \caption{Self similar PIFS.}
  \label{fig=liness}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{fracss.ps}}
  \caption{Self similar PIFS.}
  \label{fig=fracss}
\end{figure}

We can greatly increase the curves that can be described by allowing
less than the full curve to be mapped to each half.  For example, if
instead of the full interval, only the first two-thirds are mapped
onto the two halves, the resulting fractals are not self-similar,
figures~\ref{fig=linens} and~\ref{fig=fracns}.
The shape of the full
fractals is not found back in smaller versions.

As another example, the ``shock-expansion'' train of
figure~\ref{fig=train} was created by mapping the two largest
oblique line segments into themselves as in figure~\ref{fig=liness}.
The outer quarters are reduced copies of the half they are in.
Hence this train required 6 mappings:  two each for the largest
two oblique line segments and one each for the outer quarters.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{linens.ps}}
  \caption{Arbitrary PIFS.}
  \label{fig=linens}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{fracns.ps}}
  \caption{Arbitrary PIFS.}
  \label{fig=fracns}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{train.ps}}
  \caption{Arbitrary PIFS.}
  \label{fig=train}
\end{figure}

In general, we can partition the domain of interest and for each piece
search for a matching larger region that is similar in shape.  The
resulting mappings create a partitioned IFS, or PIFS, (called a local
IFS in \cite{barnsley}).  PIFS are a much more powerful way of
approximating curves than self-similar IFS.  For example, smooth
curves can be broken up into pieces small enough that they can be
approximated by, say, straight line segments: these straight line
segments can then be described as reduced copies of larger straight
line segments.  If the curve is merely continuous, it can be broken up
into segments that are approximately constant.

However, PIFS representations do not only allow a description of smooth
curves that can be split into a limited number of simple pieces, but
also of fractal curves in which similar features tend to repeat.  For
turbulent flows, a segment of the flow field may resemble some larger
segment elsewhere, even though the similarity only extends over a
limited range.

\section{Smoothed iterated function systems}
\label{sec=spifs}

The partioned iterated function systems discussed so far
tended to produce jump discontinuities, while turbulent flows
are believed to be smooth for finite Reynolds number.
However, like for the Cantor set before, smoothing can be used
to prevent the scales from actually becoming infinitely small.
This will produce solutions that are actually smooth viewed
on very small scales.

A simple smoothing is to convolve the result of the PIFS each
iteration with a smoothing function.  If the smoothing function is
positive, smoothing does not affect the convergence, and the result
will be a smooth version of the fractal.

In order to reduce the effects of the smoothing on the larger scales,
a higher-order mixed-sign smoothing can be used.  The extreme example
of such smoothing is Fourier truncation, in which all wavenumbers
below some maximum value $k=K$ are retained unchanged and larger
wavenumbers eliminated.  Despite the fact that Fourier-truncation
smoothing may in theory affect the convergence of the IFS, this was
never observed in any of the examples we computed.

Figures~\ref{fig=trains1} through \ref{fig=trains8} show the effect
of Fourier truncation smoothing.   Note the Gibbs phenomenon
typical in Fourier-truncation smoothing.  This would not occur
for smoothing with a positive smoothing function.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{trains1.ps}}
  \caption{The fractal of figure~\ref{fig=train} smoothed using
    Fourier truncation at wavenumber $k=2$.}
  \label{fig=trains1}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{trains2.ps}}
  \caption{The fractal of figure~\ref{fig=train} smoothed using
    Fourier truncation at wavenumber $k=4$.}
  \label{fig=trains2}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{trains4.ps}}
  \caption{The fractal of figure~\ref{fig=train} smoothed using
    Fourier truncation at wavenumber $k=16$.}
  \label{fig=trains4}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{trains8.ps}}
  \caption{The fractal of figure~\ref{fig=train} smoothed using
    Fourier truncation at wavenumber $k=256$.}
  \label{fig=trains8}
\end{figure}

\chapter{Fractal Representation of Velocity Fields}
\label{sec=1d}

The next step to be addressed is the construction of smoothed,
partitioned iterated function systems (SPIFS) for turbulent velocity
fields.  The basic procedure can be summarized as:
\begin{itemize}
\item Divide the flow into regions;
\item For each region, try to find a matching larger region
with a similar velocity field, save for a change of scale and
and an additional translational velocity;
\item If no match can be found for a region, divide it into
smaller pieces.
\end{itemize}
More succinctly, the key is to find matching larger and smaller
flow regions.  We will loosely refer to those flow regions
as eddies, the smaller region being the target eddy, and the larger
region that resembles it the source eddy.

According to the contraction mapping principle, up to about the error
that the larger regions approximate the smaller regions, the found IFS
will reproduce the given flow field.  More interesting is what happens
on scales smaller than those forced to be correct.  By construction,
the IFS will also represent those as reduced copies of larger scales.
This follows from the fact that the final velocity field results from
repeated application of the IFS machine.  During each application, the
larger scale shapes are further propagated into the domain of the
smaller scales.

The hope is that such short scales will provide an approximation to a
typical possible realization of the small scales in actual flows.  The
added translational velocity would presumably not make a difference on
behalf of the Galilean invariance of the equations of motion.

The change of scale does cause a change in Reynolds number, which
will make a difference.  However, we propose to use smoothing to
simulate the Reynolds number effects on the short scales.
Unless this can be done very well, the representation at Kolmogorov
length scales should be poor.  Hence, it can be expected that
IFS representations of this type will be more effective for high Reynolds
number applications with a significant inertial range.

Note that the smoothing of the short scales is important, since it
represents the Reynolds number dependence of the turbulence in the
absence of other effects such as walls. However, it may not be
important to simulate every small detail of the viscous processes. The
most important effect to simulate could be the dissipation of kinetic
energy into heat. To do so, the IFS would provide estimates for the
small scales near the viscous regime.  On the other hand, the shape of
the larger eddies could have an important Reynolds dependence which
could require modeling too.

The partitioning makes handling of flows with regions of
widely varying turbulent properties at least in principle conceivable.
An eddy in a part of the flow resembling homogeneous turbulence is
more likely to resemble a somewhat larger homogeneous
turbulence eddy than a somewhat larger eddy in some highly sheared
part of the flow field.  Similarly eddies in flows with
directionality, proximity of the wall, etcetera, are more likely to be
matched by eddies with the same orientation, scaled distance from the
wall, ...  A consequence of this is that unresolved small scales in
parts of the flow that behave, say, like homogeneous turbulence are
constructed from larger homogeneous turbulence eddies, not from eddies
in flow regions with very different properties.

From such ideas, it becomes clear that for a given flow region and its
larger image, the more nontrivial the match between the two regions is,
the better the IFS representation is likely to be.  The current study
will restrict itself to the representation of one-dimensional velocity
traces to keep things simple.  However, it is clearly much more likely
to have a chance resemblance of one-dimensional cuts through two
eddies than of the complete eddy shapes.  Hence, the present study is
limited.

It should be a much better idea to find matching eddies using their
complete three-dimensional flowfields. Even better, four-dimensional
space-time regions could be matched, so that not only the full shape,
but also the time evolution produces a close match
between the target and source eddies.

However, manipulating a three- or four-dimensional DNS or LES data
set is a rather formidable task for an exploratory study such
as this one. Hence the much more limited scope of the present
investigation.

\section{Data to be represented fractally}
\label{sec=data}

Two sets of turbulent flow data were used in this study.  The first
set consists of decaying periodic-box turbulence.  These data were
created by a periodic-box, spectral, compressible, direct
numerical simulation (DNS) code.  The flow is essentially the same as
that computed by Erlebacher {\em et al} \cite{erlebacher}.  The
initial conditions were taken from the grid turbulence measurements of
Comte-Bellot and Corrsin \cite{comte}, but scaled up to a RMS Mach
number of 0.3 (in a very small periodic box) and then advanced in time
using a spectral collocation method.  Unlike Erlebacher {\em et al},
no adjustments were made to reduce the CFL limitations at low Mach
number; also, the Sutherland viscosity law was used.  Despite these
minor differences, our results agree with those of \cite{erlebacher}
within the uncertainty in the initial conditions.  The current results
extend over 3/4 of the range computed by Erlebacher {\em et al}
because of queue time limits.  The code ran on a Thinking Machines
CM-2 SIMD hypercube computer at a scaled speed of almost one GFlop
using a $128^3$ mesh.

The initial velocity field was initialized to random numbers and then
scaled down to the experimental spectrum of Comte-Bellot and Corrsin.
The initial random velocity has a directionality along the mesh
directions, since a uniform rather than a normal probability
distribution was used for the velocity components.  However, the
effect seems minor and nonlinear effects will presumably eliminate it
quickly for all but the slowest modes.

Despite the fact that the initial spectrum agrees well with the
experiments, figure~\ref{fig=e11}, the initial turbulence does not.
This indicates that the phase information not captured by the spectrum
is significant.  For example, the initial vorticity, dissipation, and
kinetic energy seem much less ``organized'' than those produced by the
Navier-Stokes equations at the final computed time.  Compare the
two-dimensional flow cross sections figures~\ref{fig=ome20}
through~\ref{fig=ekin}.

The probability density distribution of the initial strain rates shows
equal likelihood of positive and negative mid strain rates,
figure~\ref{fig=strain0}.  This is expected regardless of the precise
details of our initialization: in the random initial data opposite
flow fields have the same chance of occurring.  However, the results
for later times show the well known propensity of the Navier-Stokes
equations to produce a positive mid strain rate,
figure~\ref{fig=strain}.  (Figures~\ref{fig=strain0}
and~\ref{fig=strain} were plotted using an oblique axis system, with
the largest strain rate horizontal and the mid strain rate under 60
degrees so that the final shape of the distribution is independent of
the ordering of the strain rates and so that opposite flow fields show
up as points symmetric about the horizontal.)

Various statistics, figures~\ref{fig=re} to~\ref{fig=fl}, show similar
``transient'' effects, in which the flow adjusts from the random
initial velocity to the usual behavior that the Navier-Stokes
equations tend to produce.  This is most clear from the ``constant''
in the $k,\epsilon$-model, figure~\ref{fig=ce2}.  The tendency of the
Navier-Stokes equations to prefer certain flow properties might be
compared to an approach to a ``strange attractor'', \cite{peitgen},
though actually this flow eventually evolves to a state of rest
instead.  Hopefully, it will make fractal representation more
advantageous.  Certainly the random initial velocity of
figure~\ref{fig=ekin0} does not look very encouraging for attempts to
find similarities between different regions of significant size.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{e11.ps}}
  \caption{1-D spectra at various times}
  \label{fig=e11}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{ome20.ps}}
  \caption{Vorticity for the random initial velocity}
  \label{fig=ome20}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{ome2.ps}}
  \caption{Vorticity for the final computed velocity}
  \label{fig=ome2}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{eps0.ps}}
  \caption{Dissipation for the random initial velocity}
  \label{fig=eps0}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{eps3.ps}}
  \caption{Dissipation for the final computed velocity}
  \label{fig=eps3}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{ekin0.ps}}
  \caption{Energy for the random initial velocity}
  \label{fig=ekin0}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{ekin.ps}}
  \caption{Energy for the final computed velocity}
  \label{fig=ekin}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 350 582 760]{strain0.ps}}
  \caption{Principal strain rates for the random initial velocity}
  \label{fig=strain0}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 350 582 760]{strain.ps}}
  \caption{Principal strain rates for the final computed velocity}
  \label{fig=strain}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{re.ps}}
  \caption{Taylor microscale Reynolds number}
  \label{fig=re}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{u.ps}}
  \caption{RMS velocity}
  \label{fig=u}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{eps.ps}}
  \caption{Dissipation}
  \label{fig=eps}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{ce2.ps}}
  \caption{k-$\epsilon$ model constant}
  \label{fig=ce2}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{sk.ps}}
  \caption{Velocity derivative skewness}
  \label{fig=sk}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{fl.ps}}
  \caption{Velocity derivative flatness}
  \label{fig=fl}
\end{figure}

The other data set used in this study was kindly made available by
Alan Wray from the Stanford/NASA Ames Center of Turbulence Research.
These data consisted of forced periodic-box turbulence on a $512^3$
mesh.  Since Wray's results are for a significantly higher Taylor
microscale Reynolds number, 169 rather than about 40, with a more
significant inertial range, it was hoped that their fractal
representation would be much better.  However, results so far
are unclear.

\section{Fractal representation procedure}
\label{sec=1dproc}

\subsection{Nature of the mapping}

Figure~\ref{fig=exammap} shows an example SPIFS representation of a
one-dimensional velocity trace from the $128^3$ point DNS simulation
described above.  The top curve shows the trace
of the velocity field according to the DNS.  The lower part of the
figure shows its partitioning into target regions.  Each target region
is represented as a reduced copy of a source region, whose extend is
indicated by the horizontal ranges in the top part of the figure.  The
highest indicated target range corresponds to the first target, and
the start of the range is indicated by the vertical tick mark.
Each mapping adds a vertical shift (translational velocity) and
a reduction in the velocity variations.

Applying the shown set of mappings to the top trace results in the
bottom trace of figure~\ref{fig=exammap}. It is a reasonable,
though not highly accurate match to the original.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{exammap.ps}}
  \caption{An example partitioned mapping.}
  \label{fig=exammap}
\end{figure}

%\begin{figure}[htbp]
%  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{exammapu.ps}}
%  \caption{The same mapping without the smoothing.}
%  \label{fig=exammapu}
%\end{figure}

If the IFS mapping is applied repeatedly, starting from any arbitrary
initial distribution, it converges to the bottom trace of
figure~\ref{fig=examcmp}.  Since the original curve is not needed to
generate the bottom curve, all information needed for it is in the
mappings. In other words, we have approximated our initial curve
by a set of mappings alone.

In addition to the restriction to one-dimension and Reynolds number
effects, the main limitations of the current implementation are the
small number of source regions actually considered and the simplistic
check for similarity between the source and target regions.  These
limitations will be discussed subsequently.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{examcmp.ps}}
  \caption{The final SPIFS.}
  \label{fig=examcmp}
\end{figure}

\subsection{Selection of the targets and possible sources}

To create the target segments in figure~\ref{fig=exammap}, the curve
is first split into two halves.  For each half, a matching larger
segment of curve is sought to act as a source.  The term ``larger'' is
here understood to mean that the velocity variation of the source
segment is reduced to produce the target. It is allowable that the
spatial extend of the source is actually less than that of the
target.

If no suitable source can be found for a target, it is subdivided into
two and a source is sought for each half.  It follows that the
possible target regions are limited; for example, the length of any
target segment is always a negative power of two.  However, this is
not a very important restriction: if an arbitrary region can be
described satisfactorily by a certain map, so can the regions found by
the binary division process that are in it.

More importantly, in the present implementation the algorithm only
takes a very limited set of source regions into account.  In
particular, the length of the source segment is presently restricted
to be an integer multiple of the length of the target, or vice-versa.
Thus only eddies which differ by an integer multiple in size are
considered.

The reason for this limitation is that the present procedure performs
an exhaustive search for source regions.  The more appropriate
search procedure would probably be simulated annealing and/or
a genetic algorithm.  This will be added at a later time.

\subsection{Identification of the best sources}

After potential source and target regions have been determined, they
need to be compared in order to identify the best match. To do so,
first the additional translational velocity in the map and the scaling
factor of the velocity variations must be determined.  This is done by
a least square minimization of the differences between target and
mapped source.

The minimization is subject to the constraint that the
map must reduce the variations in velocity. Typically, at least a 15\%
reduction in velocity variations is demanded.  In case the horizontal
extend of the source is less than or equal to that of the target, it
was found that better SPIFS approximations result from requiring that
the velocity variations are reduced by at least 50\%.  In other words,
a significant reduction in at least one direction is desired.

To speed up the search, initially only two points are compared to see
whether a source can represent a given target to a chosen accuracy.
If the source passes that test, the number of points being compared is
successively doubled until the full number of mesh points in the
target region is reached, or the desired accuracy can no longer be
met.

If more than one source meets the demanded accuracy, the best is
selected. If none meets it, the target is divided into two and the
process is repeated.

\chapter{Results}
\label{sec=results}

\section{Approximation of the viscous effects}

As mentioned in the introduction, even for the inviscid scales of
turbulence some approximation of the viscous effects is required to
simulate Reynolds number dependence effects.

To generate the smoothed PIFS of figures~\ref{fig=exammap}
and~\ref{fig=examcmp}, the smoothing effects of the viscous terms were
approximated by simply cutting off all wavenumbers above 32.  If this
smoothing is turned off, the obtained PIFS appears as
figures~\ref{fig=exapifsm} and \ref{fig=exapifsc}. These are clearly a
poorer approximation, hence the simple representation of the viscous
effects does help improve the obtained mapping.

%\begin{figure}[htbp]
%  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{examcmpu.ps}}
%  \caption{The final SPIFS evaluated without smoothing.}
%  \label{fig=examcmpu}
%\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{exapifsm.ps}}
  \caption{The corresponding PIFS mapping.}
  \label{fig=exapifsm}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{exapifsc.ps}}
  \caption{The final PIFS.}
  \label{fig=exapifsc}
\end{figure}

More sophisticated procedures may be formulated.  For example, after
reduction of the source, the tail end of its spectrum may be adjusted
to be the same as that of the original.  So far, attempts in this
direction have not been very successful.  One difficulty is that the
spectrum only represents the magnitude of the Fourier modes, while
there is nontrivial information in the phase.  This is seen, for
example, from the DNS results in section~\ref{sec=data}, where the
initial data (random phase) differ greatly from the later times
(Navier-Stokes produced phase.)

Results so far seem to indicate that merely modifying the spectrum of
the source after reducing it tends to reduce the regions of steep
gradients, while what is really desired is a reduction of the small
scale processes.  More study is needed how Reynolds number effects
modify the shape of the eddies, and not just the amplitudes of the
modes.

The contraction mapping principle is no longer formally valid when the
Reynolds number effects are represented using a smoothing as done
here.  Yet, no lack of iterative convergence was ever observed.

\section{Desirability of expansive mappings}

In the application of PIFS to image compression, \cite{barnsley}, the
practice is to map larger picture regions to smaller ones.  However,
for velocity traces it was found advantageous also to allow
smaller regions to be mapped onto larger ones. The reason is
that the velocity traces tend to show an ``intermittency'' where
quite extensive smooth regions alternate with shorter regions of rapid
variation.  An effective way to create smooth regions is to stretch
short regions; the stretching smooths the velocity distribution.

Stretching source regions will not affect convergence of the IFS as
long as the velocity variations of the target remain less
than those of the source. Instead, increasing the number of
sources considered can only reduce the deviations between target
and mapped source, which reduces the maximum possible error in
the final SPIFS correspondingly.

To verify that allowing stretched source regions does improve the
representation, stretching was turned off in
figures~\ref{fig=exammapr} and \ref{fig=examcmpr}.  Although the number
of required mappings did not change significantly, clearly the
resulting representation is much worse than in
figures~\ref{fig=exammap} and \ref{fig=examcmp}.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{exammapr.ps}}
  \caption{The SPIFS mapping without expansions.}
  \label{fig=exammapr}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{examcmpr.ps}}
  \caption{The SPIFS without expansions.}
  \label{fig=examcmpr}
\end{figure}

\section{Accuracy of the representation}
\label{sec=accuracy}

If greater accuracy is desired, it can be demanded that the source
segments approximate the target segments more closely.  This
requires that the targets are further subdivided, so that the number of
mappings increases.  It can be show from the contraction mapping
principle that any accuracy can be achieved for sufficiently many
mappings.

Figure~\ref{fig=examres} shows the convergence of the SPIFS-generated
trace (broken line) to the actual velocity trace (solid line) when the
accuracy requirements are increased.  Table~\ref{tab=examres} shows
the corresponding error in velocity, average and RMS velocity, as well
as the RMS velocity gradient and its skewness and flatness.  It
verifies the theoretically expected convergence.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{examres1.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{examres2.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{examres3.ps}}
  \caption{The effect of increasing accuracy requirements.}
  \label{fig=examres}
\end{figure}

Table~\ref{tab=examres2} shows the same for the trace half a period
distant in the two normal directions.

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
maps & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      14 & 0.077 & $-$19.0 & 55.1 & 0.0630 & 351. & $-$0.16 & 7.82 \\
      16 & 0.043 & $-$19.0 & 55.5 & 0.0621 & 351. & $-$0.24 & 7.87 \\
      28 & 0.016 & $-$19.0 & 55.7 & 0.0619 & 358. & $-$0.05 & 7.15 \\
      35 & 0.008 & $-$19.0 & 55.7 & 0.0619 & 360. & $-$0.06 & 6.96 \\
$\infty$ &   0   & $-$19.0 & 55.8 & 0.0619 & 362. & $-$0.10 & 6.97 \\
\hline
\end{tabular}
\end{center}
\caption{The effects of increasing accuracy requirements; parallel
velocity component, decaying turbulence, first trace of mesh.}
\label{tab=examres}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
maps & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
       9 & 0.102 & $-$14.9 & 65.9 & 0.0310 & 249. & 1.54 & 8.83 \\
      16 & 0.049 & $-$14.7 & 67.3 & 0.0294 & 275. & 0.76 & 5.62 \\
      22 & 0.019 & $-$14.8 & 67.4 & 0.0294 & 279. & 0.59 & 5.04 \\
      34 & 0.007 & $-$14.8 & 67.5 & 0.0295 & 280. & 0.59 & 4.70 \\
$\infty$ &   0   & $-$14.8 & 67.6 & 0.0294 & 281. & 0.60 & 4.64 \\
\hline
\end{tabular}
\end{center}
\caption{The effects of increasing accuracy requirements; parallel
velocity component, decaying turbulence, middle trace of mesh.}
\label{tab=examres2}
\end{table}

Figures~\ref{fig=wrayres} and \ref{fig=wrayres2}, and tables
\ref{tab=wrayres} and \ref{tab=wrayres2}, show the convergence
for the forced turbulence of Wray, for the velocity components
parallel and normal to the trace.  Again the procedure converges.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{wrayres1.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{wrayres2.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{wrayres3.ps}}
  \caption{The effect of increasing accuracy requirements;
    parallel velocity component, data of Wray.}
  \label{fig=wrayres}
\end{figure}


\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
maps & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      11 & 0.090 & $-$2.86 & 3.37 & 0.361 & 16.2 & 0.19 & 14.1 \\
      22 & 0.021 & $-$2.86 & 3.56 & 0.322 & 20.6 & 0.43 & 3.76 \\
      34 & 0.009 & $-$2.86 & 3.57 & 0.321 & 21.2 & 0.35 & 3.49 \\
      52 & 0.005 & $-$2.86 & 3.57 & 0.321 & 21.3 & 0.40 & 3.53 \\
$\infty$ &   0   & $-$2.86 & 3.57 & 0.321 & 21.4 & 0.42 & 3.64 \\
\hline
\end{tabular}
\end{center}
\caption{The effects of increasing accuracy requirements;
    parallel velocity component, data of Wray.}
\label{tab=wrayres}
\end{table}


\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{wrayres5.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{wrayres6.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{wrayres7.ps}}
  \caption{The effect of increasing accuracy requirements;
    normal velocity component, data of Wray.}
  \label{fig=wrayres2}
\end{figure}


\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
maps & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      13 & 0.061 & $-$0.402 & 5.33 & 0.00421 & 29.2 & 0.25 & 7.49 \\
      26 & 0.021 & $-$0.381 & 5.34 & 0.00389 & 27.1 & 0.24 & 6.05 \\
      42 & 0.010 & $-$0.384 & 5.35 & 0.00392 & 28.1 & 0.20 & 4.94 \\
      48 & 0.008 & $-$0.385 & 5.35 & 0.00393 & 28.1 & 0.22 & 4.79 \\
$\infty$ &   0   & $-$0.386 & 5.36 & 0.00395 & 28.4 & 0.19 & 4.70 \\
\hline
\end{tabular}
\end{center}
\caption{The effects of increasing accuracy requirements;
    normal velocity component, data of Wray.}
\label{tab=wrayres2}
\end{table}


The number of mappings required to represent the data may be reduced
by simultaneously mapping several velocity traces.  This is demonstrated
in table~\ref{tab=multig} for the velocity component parallel to the
traces and table~\ref{tab=multigv} for the normal component.  For example,
in the top of table~\ref{tab=multig}, the total number of maps to represent
two traces to a chosen accuracy is 75, of which 34 are to represent
the first trace.  But if for each trace also source segments taken from
the other trace are allowed, the number of mappings reduces to 71,
and 32 for the first trace.

Mapping a total of 16 traces simultaneously, the number of mappings
needed for the first trace reduces to 23.  The traces used in
tables~\ref{tab=multig} and \ref{tab=multigv} were chosen uniformly
spaced in a single flow cross section.  The results indicate that
significant reductions in the number of mappings can be achieved
by extending the data base.

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{lcccccccc}
\hline\hline
method & maps & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ &
$\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
\multicolumn{9}{c}{2 traces} \\
\hline
single &   34/75  & 0.011 & $-$3.33 & 5.58 & 0.182 & 24.5 & 0.41 & 4.38 \\
global &   32/71  & 0.012 & $-$3.33 & 5.58 & 0.183 & 24.6 & 0.39 & 4.34 \\
exact  & $\infty$ &   0   & $-$3.33 & 5.59 & 0.182 & 24.9 & 0.39 & 4.08 \\
\hline
\multicolumn{9}{c}{4 traces} \\
\hline
single &  34/152  & 0.012 & $-$3.25 & 5.53 & 0.197 & 27.9 & 0.92 & 6.43 \\
global &  28/136  & 0.011 & $-$3.25 & 5.53 & 0.197 & 27.8 & 0.92 & 6.34 \\
exact  & $\infty$ &   0   & $-$3.25 & 5.53 & 0.197 & 28.2 & 1.03 & 7.75 \\
\hline
\multicolumn{9}{c}{8 traces} \\
\hline
single &  34/312  & 0.012 & $-$3.34 & 5.40 & 0.217 & 26.7 & 0.75 & 5.82 \\
global &  25/253  & 0.015 & $-$3.34 & 5.40 & 0.217 & 26.7 & 0.73 & 5.77 \\
exact  & $\infty$ &   0   & $-$3.34 & 5.41 & 0.217 & 27.1 & 0.84 & 6.79 \\
\hline
\multicolumn{9}{c}{16 traces} \\
\hline
single &  34/616  & 0.016 & $-$3.14 & 5.54 & 0.193 & 26.4 & 0.50 & 5.82 \\
global &  23/470  & 0.015 & $-$3.14 & 5.54 & 0.193 & 26.4 & 0.49 & 5.70 \\
exact  & $\infty$ &   0   & $-$3.14 & 5.54 & 0.193 & 26.9 & 0.58 & 6.55 \\
\hline
\end{tabular}
\end{center}
\caption{Group compression of traces,
    parallel velocity component, data of Wray.}
\label{tab=multig}
\end{table}


\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{lcccccccc}
\hline\hline
method & maps & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ &
$\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
\multicolumn{9}{c}{2 traces} \\
\hline
single &   42/95  & 0.012 & 1.69 & 5.27 & 0.130 & 34.4 & 0.72 & 9.34 \\
global &   37/71  & 0.013 & 1.70 & 5.27 & 0.130 & 34.5 & 0.70 & 9.33 \\
exact  & $\infty$ &   0   & 1.69 & 5.27 & 0.130 & 35.0 & 0.72 & 9.57 \\
\hline
\multicolumn{9}{c}{4 traces} \\
\hline
single &  42/193  & 0.017 & 0.123 & 5.69 & 0.0900 & 38.2 & 0.58 & 8.20 \\
global &  33/166  & 0.017 & 0.123 & 5.69 & 0.0900 & 38.4 & 0.57 & 8.24 \\
exact  & $\infty$ &   0   & 0.123 & 5.69 & 0.0899 & 39.3 & 0.61 & 9.63 \\
\hline
\multicolumn{9}{c}{8 traces} \\
\hline
single &  42/377  & 0.017 & 0.288 & 5.31 & 0.0980 & 36.3 & 0.34 & 7.03 \\
global &  31/308  & 0.018 & 0.289 & 5.31 & 0.0980 & 36.4 & 0.31 & 6.94 \\
exact  & $\infty$ &   0   & 0.288 & 5.31 & 0.0979 & 37.3 & 0.33 & 8.03 \\
\hline
\multicolumn{9}{c}{16 traces} \\
\hline
single &  42/749  & 0.025 & 0.0834 & 5.30 & 0.0900 & 36.9 & 0.15 & 7.08 \\
global &  28/571  & 0.023 & 0.0832 & 5.30 & 0.0900 & 36.8 & 0.14 & 7.12 \\
exact  & $\infty$ &   0   & 0.0835 & 5.30 & 0.0899 & 37.8 & 0.15 & 8.08 \\
\hline
\end{tabular}
\end{center}
\caption{Group compression of traces,
    normal velocity component, data of Wray.}
\label{tab=multigv}
\end{table}


\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{multi2g.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{multi4g.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{multigg.ps}}
  \caption{Group compression of 2, 4, and 16 traces,
    parallel velocity component, data of Wray.}
  \label{fig=multig}
\end{figure}


\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{multi2gv.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{multi4gv.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[30 35 682 370]{multiggv.ps}}
  \caption{Group compression of 2, 4, and 16 traces,
    normal velocity component, data of Wray.}
  \label{fig=multigv}
\end{figure}


\section{Data compression}
\label{sec=compress}

The question to ask at this time is how efficiently the data are
represented by the fractal representation compared to a conventional
representation.  After all, one goal is to obtain a
valid representation of high Reybolds number turbulence
using a manageable number of degrees of freedom.

Even for the lower Reynolds numbers that allow direct numerical simulation,
such as those studied here, it would
be desirable to obtain a reduced description.
A typical DNS data set tends to be
big, and an effective compression algorithm would certainly help make
it easier to archive or transmit the results.
It may be noted that PIFS have been quite successful doing
this for image compression\cite{barnsley}.

The conventional representation that most accurately describes
periodic-box data is in terms of Fourier modes.  This assumes that the
errors are measured in least square sense, and follows from a Proper
Orthogonal Decomposition argument\cite{holmes}.  Hence we will compare
the SPIFS results with Fourier truncation representation.  We will
also compare with the more easily obtained but suboptimal Fourier
interpolate of sampled data.  We will restrict ourselves to
some samples,
leaving an exhaustive comparison to a later time at which the
best procedures are better understood.

First, as shown in the previous section, a SPIFS with a relatively
small number of maps can reproduce a velocity trace quite accurately.
As a reference, figure~\ref{fig=examint} shows the 32 point
Fourier interpolate for the same data as figure~\ref{fig=examcmp}.
It is not as good as the 14 map representation.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{examint.ps}}
  \caption{The 32 point Fourier interpolate for the same data
    as figure~\ref{fig=examres}.}
  \label{fig=examint}
\end{figure}

Unfortunately, a data point only requires one floating point
number, while each map requires at least two: the translational velocity
and the reduction factor. Further the maps also involve some
integers, being the subdivision level of the target and its
position, the reduction factor, and the starting point of the source.
The first three are small integers, but the starting point of the
source should probably really be counted as equivalent to
a real number.  If we rather arbitrarily put one mapping equivalent
to four mesh point values, a 32 point interpolate would
correspond to only 8 maps.  Also, the best representation
is not obtained by the Fourier interpolate, but by the Fourier
truncate.

Taking four points (or rather, modes) for each mapping,
the Fourier truncation results corresponding to
table~\ref{tab=examres} are shown in table~\ref{tab=examfout}.
The second column shows that the maximum error in the solution is
in each case smaller for the Fourier truncate than for the SPIFS.
Yet it is somewhat encouraging that the difference are not that large
until the number of modes increases and the exponential
convergence of the Fourier truncate takes over.  After all, the
main intended application is for high Reynolds number flows where
the number of modes would be far too small to reach the exponential
range.

Further, in the high Reynolds number case the hope is not so much for
exact reproduction of a particular instance of the small scales, but
of an instance with typical small-scale statistics.  The low resultion
SPIFS case compares favorably with the Fourier truncate on the
statistics of the velocity gradient.  Similarly, the SPIFS compares
well with the Fourier interpolate, table~\ref{tab=examfoua}.

The SPIFS seems noticably better than the linear interpolate using a
corresponding number of points, table~\ref{tab=examlini}.  Note that
the average velocity and average square velocity are not
representative for the errors in the velocity for the linear
interpolation case because of an effect of the periodicity;
representative errors in velocity are found in the $\delta u$ column.

However, this particular velocity trace has a thin region
with a large velocity gradient, as seen in figure~\ref{fig=examint}.
This trace does not seem to be typical.  Another trace, half a period
away, compares unfavorably with the Fourier truncate, as shown
in tables~\ref{tab=examres2} and ~\ref{tab=examfou2t}.  It does
perform better than linear interpolation, table~\ref{tab=examlin2}.

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      55 & 0.055 & $-$19.0 & 55.6 & 0.0622 & 344. & $+$0.07 & 5.96 \\
      63 & 0.027 & $-$19.0 & 55.7 & 0.0620 & 357. & $+$0.24 & 6.47 \\
     111 & 0.001 & $-$19.0 & 55.8 & 0.0619 & 362. & $-$0.11 & 7.03 \\
$\ge$128 &   0   & $-$19.0 & 55.8 & 0.0619 & 362. & $-$0.10 & 6.97 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in Fourier truncation for the same data as
table~\ref{tab=examres}}
\label{tab=examfout}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      32 & 0.404 & $-$19.4 & 54.8 & 0.0666 & 275. & $+$1.01 & 3.94 \\
      64 & 0.050 & $-$19.0 & 55.8 & 0.0616 & 362. & $+$0.21 & 6.97 \\
$\ge$128 &   0   & $-$19.0 & 55.8 & 0.0619 & 362. & $-$0.10 & 6.97 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in Fourier interpolation for the same data as
table~\ref{tab=examres}}
\label{tab=examfoua}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      56 & 0.173 & $-$19.1 & 55.8 & 0.0667 & 319. & $-$0.06 & 5.31 \\
      64 & 0.147 & $-$19.0 & 55.8 & 0.0653 & 326. & $-$0.12 & 6.31 \\
     112 & 0.044 & $-$19.0 & 55.8 & 0.0624 & 348. & $-$0.01 & 7.01 \\
$\infty$ &   0   & $-$19.0 & 55.8 & 0.0619 & 362. & $-$0.10 & 6.97 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in linear interpolation for the same data as
table~\ref{tab=examres}}
\label{tab=examlini}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      35 & 0.039 & $-$14.8 & 67.4 & 0.0295 & 262. & 0.89 & 4.26 \\
      63 & 0.017 & $-$14.8 & 67.5 & 0.0294 & 276. & 0.59 & 4.35 \\
      87 & 0.003 & $-$14.8 & 67.6 & 0.0294 & 281. & 0.59 & 4.63 \\
$\ge$128 &   0   & $-$14.8 & 67.6 & 0.0294 & 281. & 0.60 & 4.64 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in Fourier truncation for the same data as
table~\ref{tab=examres2}}
\label{tab=examfou2t}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      32 & 0.095 & $-$15.5 & 67.4 & 0.0318 & 257. & 0.90 & 4.02 \\
      64 & 0.029 & $-$14.8 & 67.6 & 0.0293 & 279. & 0.66 & 4.38 \\
$\ge$128 &   0   & $-$14.8 & 67.6 & 0.0294 & 281. & 0.60 & 4.64 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in Fourier interpolation for the same data as
table~\ref{tab=examres2}}
\label{tab=examfou2a}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      36 & 0.122 & $-$14.7 & 68.0 & 0.0426 & 242. & 0.46 & 4.27 \\
      64 & 0.046 & $-$14.8 & 67.6 & 0.0347 & 263. & 0.75 & 4.69 \\
      88 & 0.024 & $-$14.8 & 67.6 & 0.0319 & 271. & 0.63 & 4.46 \\
$\infty$ &   0   & $-$14.8 & 67.6 & 0.0294 & 281. & 0.60 & 4.64 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in linear interpolation for the same data as
table~\ref{tab=examres2}}
\label{tab=examlin2}
\end{table}

For the forced DNS results of Wray, with a somewhat higher
Reynolds number, we obtain the results of tables~\ref{tab=wrayfou1t},
\ref{tab=wrayfou1a}, and \ref{tab=wraylin1} versus \ref{tab=wrayres}
and tables~\ref{tab=wrayfou2t} \ref{tab=wrayfou2a}, and  \ref{tab=wraylin2}
versus \ref{tab=wrayres2}.

From the examples, it appears that the present SPIFS procedures can
give fairly good mid-accuracy representations of turbulent flows
compared with Fourier methods, without using periodicity in an
essential way.  But the current procedures must be improved to be
useful for data compression of DNS results.  And any superiority for
high Reynolds number description is unproven.  Most likely, such
results would require significant refinement of the crude procedures
used in this study.

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      43 & 0.038 & $-$2.86 & 3.55 & 0.325 & 16.7 & 0.08 & 2.55 \\
      87 & 0.018 & $-$2.86 & 3.57 & 0.321 & 20.3 & 0.22 & 3.16 \\
     135 & 0.007 & $-$2.86 & 3.57 & 0.321 & 21.2 & 0.30 & 3.35 \\
     207 & 0.001 & $-$2.86 & 3.57 & 0.321 & 21.4 & 0.40 & 3.62 \\
$\infty$ &   0   & $-$2.86 & 3.57 & 0.321 & 21.4 & 0.42 & 3.64 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in Fourier truncation of the same data as
table~\ref{tab=wrayres}}
\label{tab=wrayfou1t}
\end{table}


\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      32 & 0.108 & $-$2.81 & 3.40 & 0.341 & 12.5 & 0.18 & 2.74 \\
      64 & 0.059 & $-$2.85 & 3.55 & 0.321 & 19.7 & 0.15 & 2.78 \\
     128 & 0.011 & $-$2.86 & 3.57 & 0.321 & 21.3 & 0.31 & 3.26 \\
     256 & 0.001 & $-$2.86 & 3.57 & 0.321 & 21.4 & 0.42 & 3.62 \\
$\infty$ &   0   & $-$2.86 & 3.57 & 0.321 & 21.4 & 0.42 & 3.64 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in Fourier interpolation of the same data as
table~\ref{tab=wrayres}}
\label{tab=wrayfou1a}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      44 & 0.087 & $-$2.90 & 3.58 & 0.331 & 14.8 & 0.58 & 3.41 \\
      88 & 0.037 & $-$2.85 & 3.56 & 0.322 & 18.6 & 0.37 & 3.24 \\
     136 & 0.020 & $-$2.86 & 3.57 & 0.322 & 20.2 & 0.40 & 3.40 \\
     208 & 0.010 & $-$2.86 & 3.57 & 0.321 & 20.8 & 0.38 & 3.41 \\
$\infty$ &   0   & $-$2.86 & 3.57 & 0.321 & 21.4 & 0.42 & 3.64 \\
\hline
\end{tabular}
\end{center}
\caption{The errors in linear interpolation of the same data as
table~\ref{tab=wrayres}}
\label{tab=wraylin1}
\end{table}


\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      51 & 0.060 & $-$0.386 & 5.33 & 0.00397 & 22.1 & 0.58 & 4.04 \\
     103 & 0.014 & $-$0.386 & 5.36 & 0.00395 & 27.4 & 0.26 & 4.64 \\
     167 & 0.005 & $-$0.386 & 5.36 & 0.00395 & 28.2 & 0.19 & 4.67 \\
     191 & 0.004 & $-$0.386 & 5.36 & 0.00395 & 28.3 & 0.18 & 4.69 \\
$\infty$ &   0   & $-$0.386 & 5.36 & 0.00395 & 28.4 & 0.19 & 4.70 \\
\hline
\end{tabular}
\end{center}
\caption{The effects in Fourier truncation for the same data as
table~\ref{tab=wrayres2}}
\label{tab=wrayfou2t}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      32 & 0.124 & $-$0.292 & 5.25 & 0.00283 & 20.5 & 0.27 & 3.47 \\
      64 & 0.069 & $-$0.371 & 5.33 & 0.00376 & 26.5 & 0.23 & 3.94 \\
     128 & 0.014 & $-$0.389 & 5.36 & 0.00398 & 28.2 & 0.27 & 4.80 \\
     256 & 0.003 & $-$0.386 & 5.36 & 0.00395 & 28.4 & 0.19 & 4.71 \\
$\infty$ &   0   & $-$0.386 & 5.36 & 0.00395 & 28.4 & 0.19 & 4.70 \\
\hline
\end{tabular}
\end{center}
\caption{The effects in Fourier interpolation for the same data as
table~\ref{tab=wrayres2}}
\label{tab=wrayfou2a}
\end{table}

\begin{table}[htbp]
\begin{center}
\leavevmode
\begin{tabular}{cccccccc}
\hline\hline
points & $\delta u$ & $\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
      52 & 0.088 & $-$0.375 & 5.25 & 0.01600 & 20.7 & 0.47 & 3.10 \\
     104 & 0.034 & $-$0.386 & 5.36 & 0.00923 & 25.3 & 0.10 & 4.51 \\
     168 & 0.018 & $-$0.387 & 5.36 & 0.00671 & 27.0 & 0.25 & 4.80 \\
     192 & 0.014 & $-$0.386 & 5.36 & 0.00619 & 27.3 & 0.21 & 4.47 \\
$\infty$ &   0   & $-$0.386 & 5.36 & 0.00395 & 28.4 & 0.19 & 4.70 \\
\hline
\end{tabular}
\end{center}
\caption{The effects in Fourier interpolation for the same data as
table~\ref{tab=wrayres2}}
\label{tab=wraylin2}
\end{table}


If the stored flow field is intended to be used to restart the DNS,
the question further arises how important it is to store the higher
modes.  The Navier-Stokes equations seem to do a good job of restoring
them if they are missing, at least at the current relatively low
Reynolds number.  To demonstrate this, a fairly severe truncation was
performed in the middle of the computation of figures~\ref{fig=e11}
to~\ref{fig=fl}.  All wave numbers above 16 were set to zero.
Although initially there is, of course, a significant effect, at later
times there is a remarkable recovery; figures~\ref{fig=e11c}
to~\ref{fig=flc}.  At least this flow field does not appear to depend
very sensitively on the initial data, for the considered range of
times.

\vspace{.4truein}
\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{e11c.ps}}
  \caption{1-D spectra at various times}
  \label{fig=e11c}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{ome2c.ps}}
  \caption{Vorticity for the final computed velocity}
  \label{fig=ome2c}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{eps3c.ps}}
  \caption{Dissipation for the final computed velocity}
  \label{fig=eps3c}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 150 582 560]{ekinc.ps}}
  \caption{Energy for the final computed velocity}
  \label{fig=ekinc}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[30 350 582 760]{strainc.ps}}
  \caption{Principal strain rates for the final computed velocity}
  \label{fig=strainc}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{rec.ps}}
  \caption{Taylor microscale Reynolds number}
  \label{fig=rec}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{uc.ps}}
  \caption{RMS velocity}
  \label{fig=uc}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{epsc.ps}}
  \caption{Dissipation}
  \label{fig=epsc}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{ce2c.ps}}
  \caption{k-$\epsilon$ model constant}
  \label{fig=ce2c}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{skc.ps}}
  \caption{Velocity derivative skewness}
  \label{fig=skc}
\end{figure}

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 365]{flc.ps}}
  \caption{Velocity derivative flatness}
  \label{fig=flc}
\end{figure}

\section{Predictive properties}

The most important question is how well a SPIFS can predict properties
of the small scales in a turbulent flow given only limited data.  Such
predictions open up the possibility of using a SPIFS to adjust the
subgrid scales model in a LES.  Based on representation of smaller
scales as reduced copies of larger scales in similarly behaving
regions, the SPIFS could lead to predictions of where there is
relatively little or much subgrid turbulence, its directional
properties, etcetera.  These predictions could then be used to adjust
the subgrid scale model, making the model constant dynamic, adding
anisotropic viscosity tensor effects, etcetera.

As a first orientation into this problem, the 128 velocity
values were reduced to a 32 point subset, and the SPIFS
algorithm applied on those.  The question was how well the
SPIFS could reproduce the properties of the full velocity
trace using only the information of these 32 points.

Note that the SPIFS can always produce predictions that are at least
as good as the Fourier interpolate: simply use the Fourier interpolate
to create more points and then represent those accurately using a
SPIFS.  The question is whether the SPIFS can do {\it better} than
that.  To do so, the SPIFS should presumably be based only on the
correct velocity values, not the interpolated ones in between.  With
only 32 points left, this obviously does not leave much points per map.

Figure~\ref{fig=extrap} shows results obtained this way for four
different velocity traces equally spaced throughout the flow field.
The broken curves are the original velocity traces, the points are the
32 velocities used to find the SPIFS, and the solid curves are the
results of the SPIFS.  Although no actual agreement in the details can
be expected or is needed, the hope is that the SPIFS result gives some
equally likely alternate realization of the flow.  However, a visual
inspection of the results suggests that there may be systematic
differences; the SPIFS appear to have more short scale oscillations
than the original traces.

The statistical properties of the traces also show noticeable
differences, as shown in table~\ref{tab=stats}.  No clear
advantage of the SPIFS above the Fourier interpolate is evident.
Clearly, further research, and maybe a higher Reynolds number,
is needed if we want to get reliable predictions.
Yet it should be kept in mind that the present algorithm
is very crude.  In the discussion section, we will indicate
some obvious significant improvements that can be made.

\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{extrap0.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{extrap1.ps}}
  \caption{SPIFSs found using a 32 point subset.}
  \label{fig=extrap}
\end{figure}
\begin{figure}[htbp]
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{extrap2.ps}}
  \centerline{\epsfxsize=5.5in\epsfbox[50 35 662 300]{extrap3.ps}}
  \addtocounter{figure}{-1}
  \caption{Continued.}
\end{figure}

\begin{table}[htbp]
\begin{center}
\leavevmode

\begin{tabular}{cccccc}
\hline\hline
$\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
 $-$19.0 & 55.8 & 0.0619 & 362. & $-$0.10 & 6.97 \\
 $+$25.7 & 71.9 & 0.0691 & 354. & $+$0.73 & 3.76 \\
 $+$21.7 & 66.2 & 0.0597 & 365. & $+$0.16 & 4.05 \\
 $-$19.5 & 74.1 & 0.0383 & 441. & $-$1.56 & 7.59 \\
\hline
\end{tabular}

\vspace{.1truein}

\begin{tabular}{cccccc}
\hline\hline
$\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
 $-$19.0 & 54.4 & 0.0649 & 271. & $+$0.87 & 3.37 \\
 $+$25.7 & 71.8 & 0.0691 & 340. & $+$0.63 & 3.61 \\
 $+$21.7 & 65.5 & 0.0609 & 320. & $+$0.35 & 2.86 \\
 $-$19.5 & 72.8 & 0.0396 & 382. & $-$0.58 & 4.10 \\
\hline
\end{tabular}

\vspace{.1truein}

\begin{tabular}{cccccc}
\hline\hline
$\overline{u}$ & $\sqrt{\overline{u^2}}$ & $\ell$ &
$\sqrt{\overline{u^2_x}}$ & $\overline{u_x^3}/\sqrt{\overline{u_x^2}}^3$ &
$\overline{u_x^4}/\sqrt{\overline{u_x^2}}^4$ \\
\hline
 $-$20.6 & 51.0 & 0.0853 & 270. & $+$0.63 & 3.58 \\
 $+$25.5 & 65.2 & 0.0818 & 382. & $+$0.42 & 4.47 \\
 $+$17.7 & 57.4 & 0.0532 & 378. & $-$0.18 & 3.89 \\
 $-$16.2 & 68.4 & 0.0313 & 482. & $-$0.88 & 8.00 \\
\hline
\end{tabular}

\caption{Statistics at four different positions using a 32 point
subset. Top is exact, middle is the Fourier interpolate prediction,
and bottom is the SPIFS prediction.}
\end{center}
\label{tab=stats}
\end{table}

\chapter{Discussion}
\label{sec=discuss}

The theoretical properties of SPIFSs and the computed examples in this
work do seem encouraging. Yet clearly significant improvements
will be needed to reach our main objective: reliable prediction
of subscale velocity fields based on the properties of the
large eddies.

To summarize the results: we are able to give accurate SPIFS
descriptions of the one-dimensional velocity traces we tested.  Yet we
could mostly have done better by using simple Fourier interpolation of
the mesh points.

It may be conjectured that the Reynolds number of the used data is
too low to expect extrapolations based on similarity to be
very effective.  But if the techniques are to be validated
in the foreseeable future, they must be made to work at moderate
Reynolds number. Further, preliminary results for the higher
Reynolds number flow field of Wray so far do not seem to show
a clear improvement.

Since the present work has a number of serious shortcomings,
it can be hoped that results can be improved significantly by
removing these shortcomings:

\begin{enumerate}

\item The current work selects individual velocity traces and tries to
represent them by a SPIFS.  A much better approximation should be
possible by mapping all velocity traces at once.  For each target
eddy, this allows a search for a source eddy through all traces, not
just one.  For the present $128^3$ mesh, that would increase the data
base of source eddies by a factor $3\ 128^2$ or $6\ 128^2$, ensuring
better, and presumably much better, matches between targets and
sources.

\item The current work only looks for source eddies with a size that
is an integer multiple (or ratio) of the target eddies.  Clearly that
is a significant restriction on the source eddies that can qualify.

\item The current work is scalar.  Clearly there is a much greater
possibility of a chance resemblance of a single velocity component
between two eddies than of all three components.  This would
presumably affect the predictive properties of the SPIFS more
than the approximation properties.

\item The current work ignores the fact that velocity traces with
negative $u_x$ are fundamentally different from those with positive
$u_x$, freely inverting its sign in the mapping.

\item The current work is one-dimensional.  Clearly there is a much
greater possibility of a chance resemblance of two one-dimensional
cuts through an eddy than of the full three-dimensional spatial, or
four dimensional spatial-temporal velocity field.

\item The current work ignores the pressure in comparing eddies.
Comparison of the pressure field may further help determine the best
match between target and source, though at least in the
four-dimensional case the possible errors in pressure are already
restricted by the agreement in the velocity fields.

\item Since dealing with four dimensional sets could be daunting,
it would be attractive to find a one-dimensional system to test
various ideas on.  However, the one-dimensional Burgers' equation
on a periodic domain converges to a simple wave train, not resembling
turbulence, compare \cite{burgers}.  The system could possibly
be formed from the complex Ginzburg-Landau equation \cite{kuramoto}.

\item One difficulty in predicting small scale features from a limited
number of points is that the algorithm runs out of points to
compare.  Typically the difficulty is finding source eddies for target
eddies with sharp gradients.  The target eddy has not enough mesh
points to allow a division level fine enough to handle the gradients.
The current algorithm then relaxes the accuracy requirements for all
eddies.  Of course, it is not useful to globally interpolate more
points using, say, the Fourier interpolate, because the resulting
SPIFS would merely accurately reproduce that Fourier interpolate,
leaving no possibility of improvement.  But it might be acceptable to
either allow greater errors for a single failed mapping, or to use
interpolated points merely in the regions where no acceptable match
can be found.

\item The main difficulty for target eddies with steep velocity
gradients appears to be the theoretical limitation that the velocity
variations must decrease from source to target eddy.  However, it
might be possible to relax this condition under certain restrictions,
such as when it is found that the SPIFS is still contractive after
more than one application, or when the target eddy is smaller than the
source.

\item The present procedure simply models the Reynolds number effects
on the eddy shape by cutting off the tail of the spectrum.  There
should be better ways of modeling these effects.

\item The present procedure simply takes the ``best'' source eddy to
be the one with the smallest maximum difference in mesh point values
between target and mapped source. There might be a better criterion to
select the final source.  If the number of points in the target eddy
is held constant, and the size of the set of source eddies is
increased indefinitely, there will be many excellent matches of the
finite set of mesh points; which one should be chosen?

\end{enumerate}

In addition, there are many other difficulties that need to be
addressed before fractal representation will be a practical tool. We
mention:

\begin{enumerate}

\item The work involved in finding matching source eddies for each
target may be considerable.  However, due to the fact that the
velocity field is continuous, a first comparison between target and
source can start with a few points.  This is unlike fractal
compression of images \cite{barnsley}, where a few points might be
completely unrepresentative of the total.

\item This brings up the question whether the velocity field is in
fact continuous in the limit of infinite Reynolds number.

\item The work involved in evaluating the velocity field by repeated
application of the SPIFS could be prohibitive.  This evaluation
certainly does not seem trivial since it is designed to produce more
scales than we can handle in a direct computation.  Since we are
presumably mainly interested in short-scale statistics, should a
sampling of points be used, such as in the chaos game?

\item Do we need to evaluate the SPIFS? If we are mainly interested
in general features, such as the level of small scales
in a given target region, might it suffice to simply look
at the properties of the source?

\item There is at the moment no good theoretical framework to indicate
how well fractal representation can represent a turbulent flow, and
its limitations.

\item It is not clear what properties of the small scales are
the important ones to be correctly represented.

\end{enumerate}

\section*{Acknowledgments}

We thank Alan Wray for kindly making his data of forced periodic box
turbulence available, and NASA Ames and the Supercomputer Computations
Research Institute for providing various computational facilities.

\begin{thebibliography}{99}

\bibitem{barnsley}
Barnsley, M.~F.~\& Hurd, L.~P. (1992)
{\it Fractal image compression}. Peters.

\bibitem{burgers}
Burgers, J.~M. (1974)
{\it The Nonlinear Diffusion Equation: Asymptotic Solutions and
Statistical Problems}. Reidel.

\bibitem{comte}
Comte-Bellot, G. \& Corrsin, S. (1971)
``Simple Eulerian time-correlation of full and narrow-band velocity
signals in grid-generated isotropic turbulence.''
{\it J.~Fluid Mech.} {\bf 48} 273-337.

\bibitem{erlebacher}
Erlebacher, G., Hussaini, M.~Y., Speziale, C.~G., \& Zhang, T.~A. (1990)
``Towards the Large-Eddy Simulations of Compressible Turbulent Flows.''
NASA CR-187460; replaces NASA CR-178273.

\bibitem{ferziger}
Ferziger, J.~H. (1993)
``Subgrid-Scale Modeling.''
In {\it Large Eddy Simulation of Complex Engineering and Geophysical Flows},
B.~Galperin \& S.~A.~Orszag, Eds. Cambridge University Press, pp. 37-54.

\bibitem{holmes}
Holmes, P., Lumley, J.~L., \& Berkooz, G. (1996)
{\it Turbulence, Coherent Structures, Dynamical Systems and Symmetry}.
Cambridge University Press.

\bibitem{kuramoto}
Kuramoto, Y. (1997)
``Phase- and center-manifold reductions for large populations of
coupled oscillators with application to non-locally coupled
systems.''
{\it Int.~J.~Bifurcation and Chaos} {\bf 7} 789-805.

\bibitem{lesieur}
Lesieur, M. (1990)
{\it Turbulence in Fluids}.
Second revised edition, Kluwer.

\bibitem{orszag}
Orszag, S.~A., Staroselsky, I., \& Yakhot, V. (1993)
``Some basic challenges for large eddy simulation research.''
In {\it Large Eddy Simulation of Complex Engineering and Geophysical Flows},
B.~Galperin \& S.~A.~Orszag, Eds. Cambridge University Press, pp. 55-78.

\bibitem{peitgen}
Peitgen, H.-O., J\"urgen, H., \& Saupe, D. (1992)
{\it Chaos and Fractals}.  Springer-Verlag.

\bibitem{piomelli}
Piomelli, U. (1993)
``Applications of large eddy simulations in engineering: an overview.''
In {\it Large Eddy Simulation of Complex Engineering and Geophysical Flows},
B.~Galperin \& S.~A.~Orszag, Eds. Cambridge University Press, pp. 119-137.

\end{thebibliography}

\end{document}
