Wavelet Methods in Computational Fluid Dynamics ∗

This article reviews state-of-the-art adaptive, multiresolution wavelet methodologies for modeling and simulation of turbulent flows with various examples. Different numerical methods for solving the Navier-Stokes equations in adaptive wavelet bases are described. We summarize coherent vortex extraction methodologies, which utilize the efficient wavelet decomposition of turbulent flows into space-scale contributions, and present a hierarchy of wavelet-based turbulence models. Perspectives for modeling and computing industrially relevant flows are also given.

Wavelet methods in computational fluid dynamics (CFD) compose a relatively new research field and appeared in the literature only a little more than a decade ago.The ability of wavelets to identify and isolate localized structures such as shock waves and vortices, combined with the mathematical rigor of multiresolution analysis, made them attractive candidates for adaptive computational approaches and turbulence modeling.Multiresolution schemes, especially in the context of compressible Euler equations, have been intensively studied in recent years, and a number of books dealing with the subject are already available (e.g., Cohen 2000Cohen , M üller 2003)).Conversely, the use of wavelets for modeling and simulations of turbulent flows is a relatively new area of research, with a limited number of people working in the field, and the need for a systematic introduction to the subject is rather pressing.This review aims to provide an overview of wavelet methodologies in CFD, summarizes existing wavelet-based numerical algorithms, and discusses state-of-the-art turbulence modeling and simulations using wavelets, and their future potential for CFD.
CFD of turbulent flows is still a major challenge and preoccupation for the scientific community.The need to accurately predict the effect of turbulent flows impacts virtually every field of science and engineering.Turbulence is characterized by its intrinsic multiscale behavior, its self-organization into coherent structures, and a generic randomness.The major computational challenge comes from the fact that turbulence is active over a large and continuous range of length scales, which increases with Reynolds number such as Re 1/2 and Re 3/4 for two-dimensional (2D) and 3D turbulence, respectively.As a result, direct approaches that solve the governing Navier-Stokes equations at all scales without any model [i.e., by direct numerical simulation (DNS)] are limited by the available computing power, which restricts the applicability of DNS to flows of engineering interest.Hence turbulence models are unavoidable to reduce computational complexity while capturing the physics of turbulent flows.
The traditional methodology that dominated turbulence modeling over the past 40 years is the Reynolds-averaged Navier-Stokes (RANS) approach (e.g., Durbin & Reif 2001, Gatski et al. 1996), in which time-or ensemble-averaged flow fields are computed.This approach can capture some basic features in turbulent flows but cannot accurately predict the spatiotemporal characteristics.This motivated the development of large-eddy simulations (LES), in which a formal scale separation is obtained by means of low-pass filtering of the Navier-Stokes equations.The filtered Navier-Stokes equations are closed by modeling the subgrid-scale (SGS) stresses that account for ANRV400-FL42-20 ARI 9 September 2009 21:39 the effect of the unresolved small scales (e.g., see Germano et al. 1991, Lesieur & Métais 1996, Meneveau & Katz 2000).The inherent limitation of both RANS and LES approaches is that neither takes advantage of the fundamental properties of turbulence, namely, its multiscale character, its self-organization into coherent structures, and spatiotemporal intermittency.As noted by Yakhot & Sreenivasan (2005), one needs a way to detect and track energy-containing motions to describe turbulent flows with a minimum number of degrees of freedom.Wavelet multiresolution analysis offers such a capability together with a different perspective on turbulence modeling.Wavelet-based turbulence models rely on observations showing that for a given flow realization, the dynamically active scales in turbulent flows are not distributed homogeneously, neither in space nor in time, which corresponds to the flow intermittency.To benefit from this property, a suitable representation of the flow should reflect the sparsity of the fine-scale activity, in both space and time.Wavelets yield a compressed representation of turbulent flow fields; i.e., only few wavelets, with the largest in magnitude coefficients, are sufficient to represent the dynamically active part of the flow.
Wavelets offer a unique hierarchical framework for modeling and simulating turbulent flows based on the ability of wavelet multiresolution analysis to identify and isolate the energetic coherent structures that govern the dynamics of the flow.Sparse (compressed) representations of turbulent flow fields coupled with wavelet-based numerical methods allow the tight integration of the numerics and physics-based modeling.The first approach in this hierarchy is wavelet-based direct numerical simulation (WDNS), which uses wavelet-based discretizations of the Navier-Stokes equations to adapt dynamically the local resolution to intermittent flow structures.The most intriguing part of WDNS is its ability to take advantage of spatial and temporal intermittency and perform simulations of high-Reynolds number flows with computational complexity lower than commonly accepted.The next method in this hierarchy is coherent vortex simulation (CVS), which was introduced by Farge et al. (1999).The underlying idea of this method is the decomposition of the flow into coherent and incoherent contributions by means of wavelet filtering of the vorticity field.The evolution of the coherent flow is then computed deterministically, whereas the influence of the incoherent background flow is statistically modeled or neglected.To further reduce computational cost, Goldstein & Vasilyev (2004) recently proposed the stochastic coherent adaptive large-eddy simulation (SCALES).SCALES inherits CVS's ability to dynamically track the most energetic part of coherent eddies in a turbulent flow field, while (similarly to LES) modeling the effect of the less energetic (unresolved) motions.
This review is organized as follows.After a short introduction to wavelets in Section 2, we discuss different adaptive wavelet-based methods for solving the Navier-Stokes equations in Section 3, focusing on the incompressible regime.Extensions to compressible flows and techniques for applying wavelet methods to flows in complex geometries are also presented in Section 3. Section 4 provides a brief introduction to the wavelet analysis of turbulent flows with a focus on coherent vortex extraction.We review the hierarchy of wavelet-based turbulence modeling approaches in Section 5. Finally, in Section 6, we discuss the future perspectives for wavelets in CFD, in general, and turbulence modeling, in particular.

General Properties
Wavelet theory, developed mostly over the past 25 years, has generated a tremendous amount of interest in many areas of research in mathematics, physics, computer science, and engineering.The continuous wavelet transform was discovered by Grossmann & Morlet (1984)  wavelet transform by Lemarié & Meyer (1986).Daubechies (1988) developed orthogonal bases made of compactly supported wavelets, and Mallat (1989) designed the fast wavelet transform algorithm.The first application of wavelets in fluid mechanics can be traced to Farge & Rabreau (1988) (for a review, see Farge 1992).
Wavelets are basis functions that are localized in both physical space (owing to their fast decay or even compact support) and wave-number space (owing to their vanishing moments and smoothness) (see Figure 1).For comparison, the classical Fourier transform is based on functions (sines and cosines) that are well localized in wave-number space, but do not provide localization in physical space as they have global support.Due to this space/scale localization, the wavelet transform provides both spatial and scale (frequency) information, whereas the Fourier transform only provides frequency information.
A scalar field u(x) can be represented in terms of wavelet basis functions as where ū0 i and φ 0 i (x) are the scaling coefficients and functions on the lowest level of resolution, and ũμ, j i and ψ μ, j i (x) are the wavelet coefficients and basis functions of different families, μ, and levels of resolution, j.Bold subscripts denote indicies in d-dimensional space; i.e., i = (i 1 , . . ., i d ), I 0 φ , and I μ, j ψ are d-dimensional index sets, associated with scaling functions at zero level of resolution and wavelets of family μ and level j.One may think of a wavelet decomposition as a multilevel or multiresolution representation of a function, in which each level of resolution j (except the coarsest one) consists of wavelets ψ j i or a family of wavelets ψ μ, j i having the same scale but located at different positions.Scaling function coefficients, ū0 i , represent the averaged values of the field, whereas the wavelet coefficients, ũμ, j i , measure the fluctuations (details) of the field u(x) at j scale and around position of the wavelet ψ μ, j i (x).In d dimensions, there are 2 d − 1 distinctive d-dimensional wavelet families (Daubechies 1992).The above wavelet decomposition can be also applied to vector fields component-wise.
Wavelet and scaling function coefficients can be obtained, respectively, by multiplying Equation 1 by the bi-orthogonal dual wavelets ψμ, j i (x) and scaling functions φ0 i (x) and integrating over the domain, resulting in the following definition of the coefficients: where •, • denotes the L 2 -inner product, defined by u, g = u(x)g(x)d x.In the case of orthogonal wavelets, the dual wavelets are identical to primary ones, i.e., ψ = ψ.The wavelet decomposition is computationally efficient, thanks to the existence of a fast recursive transform that allows computation of wavelet coefficients from the scaling coefficients at a finer level of resolution.The fast wavelet transform yields an efficient O(N ) algorithm to compute the N wavelet coefficients ũμ, j i from the N grid point values of the function u(x) (for a more general discussion on wavelets, see Daubechies 1992, Mallat 1999).
Traditionally, 1D first-generation wavelets ψ j k are defined as translates and dilates of one basic wavelet ψ, i.e., ψ (Daubechies 1988), resulting in wavelet construction in unbounded or periodic domains with uniform sampling interval.Wavelets were further generalized to bounded domains, mainly by modifying them close to domain boundaries, but keeping them unchanged in the middle of the domain (e.g., Chiavassa & Liandrat 1997, Cohen & Daubechies 1993, Meyer 1991, Monasse & Perrier 1998).An alternative interval wavelet construction was proposed by Sweldens (1996), who constructed wavelets in physical space without the requirement of translation and dilation invariance.These wavelets, commonly referred to as second-generation wavelets, supply additional freedom to deal with arbitrary boundary conditions and irregular sampling intervals, while still being local in both space and frequency and having vanishing polynomial moments.Both first-and second-generation wavelets can be generalized easily to multiple domains as tensor product wavelets, which leads to the existence of 2 d − 1 wavelet families in d dimensions.

Wavelet Thresholding and Denoising
Filtering can be performed in wavelet space using thresholding of the coefficients.This can be considered as a nonlinear filter as the retained coefficients depend on the filtered function.The wavelet thresholding filter, [•] , is defined by where u(x) is a scalar field, > 0 stands for the nondimensional (relative) threshold parameter, and • is the norm that provides the (absolute) dimensional scaling for the filtered variable u.
For instance, in the case of velocity, the (absolute) dimensional scaling can be specified as the Once the norm is chosen, the wavelet thresholding filter (Equation 3) is uniquely defined by the nondimensional threshold parameter .For clarity of presentation, a shorter notation of the wavelet-filtered quantity u (x) = [u(x)] is used throughout the review when possible.The reconstruction error due to wavelet filtering with the nondimensional threshold parameter can be shown to be for a sufficiently smooth function u(x), where C is of order unity (Donoho 1992).This property of wavelet compression is the foundation of the active error control used in all wavelet-based numerical algorithms.Wavelet thresholding for denoising signals corrupted with Gaussian white noise is based on the fact that in wavelet space only few strong coefficients represent the signal, while the large majority of small coefficients correspond to the noise.Thus for reasonable signal-to-noise ratios, thresholding of the weak coefficients allows the recovery of a denoised signal with optimal minimal error, as shown by Donoho & Johnstone (1994).

WAVELET-BASED NUMERICAL METHODS
The mathematical properties of wavelets motivate their use for the numerical solution of PDEs.The localization of wavelets, both in scale and space, leads to effective sparse representations of functions and pseudodifferential operators (and their inverses) by performing nonlinear thresholding of the wavelet coefficients of the function and of the matrices representing the operators.The most appealing feature of wavelet analysis for the numerical solution of PDEs is the ability to estimate the local regularity of the solution, which allows self-adaptive discretizations with automatic local mesh refinement.Furthermore, the characterization of function spaces in terms of wavelet coefficients and the corresponding norm equivalences (Dahmen & Kunoth 1992, Jaffard 1992) allow diagonal preconditioning of operators in wavelet space.Finally, the existence of the fast wavelet transform yields algorithms with optimal linear complexity.For a more mathematical overview on wavelet methods for PDEs, the reader is referred to the following reviews (Cohen 2000, Dahmen 1997, M üller 2003).
Adaptation of the wavelet expansion occurs quite naturally and is based on the analysis of the wavelet coefficients of the decomposition (Equation 1).If one is interested in representing a function or a field with the fewest degrees of freedom, while still retaining a good approximation, then the wavelet thresholding criterion (Equation 3) is optimal for this purpose.However, when solving PDEs, one needs to apply an additional criterion to ensure that the wavelet basis or computational mesh is sufficient to approximate the solution throughout the time-integration step for an evolution problem or at the next iteration in the elliptic case.The adaptation strategy for an evolution problem is illustrated in Figure 2. To ensure the adequate approximation of the solution during time integration, Liandrat & Tchamitchian (1990) introduced the concept of a safety or an adjacent zone, which includes wavelets whose coefficients are, or can possibly become, significant during the period of time integration, when the grid remains unchanged.Most current wavelet-adaptation techniques are based on this strategy.In actual implementations, the safety Illustration of the adjacent zone for adaptive wavelet methods in wavelet coefficient space with scale index j and position index i.Solid pink circles indicate the positions of wavelets.The locations of wavelets with significant coefficients kept in Approximation 3 at time t are marked by the orange bell curve.The second bell-shaped region on the right of the original one marks the locations of the significant wavelet coefficients at the end of the time integration, i.e., at t + t.The safety zone is represented by red bell-shaped region.zones [neighboring wavelets at the same, one above (children) and one below (ancestors), levels of resolution] are added for each significant coefficient, and the time step is chosen to ensure that the solution does not propagate outside the safety zone.This typically results in a CFL-like constraint (Courant et al. 1928), which ensures that no energy at a given resolution scale propagates outside the safety region.The thickness of the safety zone determines the time interval during which the calculations can be carried out without modifying the computational grid.However, for computational efficiency, it was found (e.g., Schneider et al. 2006, Vasilyev 2003) that the safety zone, which includes the immediate neighboring wavelets, is the optimal.The current wavelet-based algorithms can be classified in different ways depending on whether they take full or partial advantage of wavelet analysis, namely, multiresolution properties, wavelet compression, the detection of localized structures and subsequent use for grid adaptation, fast wavelet transform, wavelet-based interpolation, and active error control.For this review, we roughly divide wavelet methods into three main classes: pure wavelet, multiresolution, and waveletoptimized methods.In addition, there are two unmistakable groups of wavelet-based methods that have unique characteristics: Lagrangian and space-time adaptive wavelet methods.Here we briefly outline the main characteristics of each class of wavelet-based algorithms.Pure wavelet schemes (either Galerkin or collocation) directly use wavelet decomposition to discretize the PDEs.Multiresolution methods were inspired by wavelets but are based on finite-difference, finite-volume, or finite-element formulation with multiresolution and automated refinement capabilities.Lagrangian wavelet methods are based either on moving wavelets or Lagrangian particle methods complemented with multiresolution capabilities.Space-time strategies include methods that either use local time stepping or treat time as an additional dimension and apply wavelet decomposition and adaptation simultaneously in space and time.Finally, wavelet-optimized methods include algorithms based on classical discretizations (e.g., by finite differences or finite volumes) but use wavelet analysis either to define adaptive meshes or to speed up the linear algebra.
In fluid dynamics, different approaches have been developed for the compressible and incompressible Euler and Navier-Stokes equations for both inert and reactive flows.The earlier attempts to solve incompressible Navier-Stokes equations using wavelet methods were based on either vorticity-stream function (Charton & Perrier 1996, Fr öhlich & Schneider 1996) or vorticity-velocity (e.g., Vasilyev & Kevlahan 2002) formulations, mainly to avoid the issues related to dynamic pressure, which limited their use to two spatial dimensions.The velocity-pressure formulation in both two and three dimensions was used by Vasilyev and colleagues (Goldstein et al. 2005, Kevlahan & Vasilyev 2005), who applied an adaptive wavelet collocation method to solve the Poisson equation to impose the incompressibility condition.Another way to impose solenoidality on the velocity field is to use divergence-free wavelets for solving the incompressible Navier-Stokes equations in two dimensions for periodic boundary conditions (Bittner & Urban 2007, Deriaz & Perrier 2006, Zhou & He 2005).In addition to model incompressible flows in complex geometries, different adaptive wavelets methods with Brinkman penalization for stationary (e.g., Schneider & Farge 2002, Vasilyev & Kevlahan 2002) and moving (e.g., Kevlahan & Vasilyev 2005) obstacles have been proposed.

Pure Wavelet Methods
This section describes pure wavelet methods of both Galerkin and collocation type.Either wavelets are used directly to discretize the underlying PDE or their properties are used for grid adaptation, preconditioning, or inversion of the operators and fast interpolation.To simplify the presentation, let us consider the general elliptic problem Lu = f , where L represents, for example, the Laplace ∇ 2 or the Helmholtz (I − ∇ 2 ) operator, which arises when solving the incompressible Navier-Stokes equations.Boundary conditions may occasionally have been incorporated into L.

Adaptive wavelet Galerkin methods.
In Galerkin methods, the unknown u is developed into a wavelet series given by Equation 3 considering only coefficients with absolute value above a given threshold (e.g., Maday et al. 1991, Maday & Ravel 1992).Substituting Approximation 3 into the differential equation and requiring that the residual vanishes with respect to test functions1 result in a system of equations for the unknown wavelet and scaling function coefficients, ũμ, j i and ū0 i , which can be solved at each time step.The underlying adaption strategy is based on analyzing the wavelet coefficients ũμ, j i and the adjacent zone as explained at the beginning of Section 3 and illustrated in Figure 2.With such a strategy approximation, Equation 3 automatically tracks the solution in space and scale.Boundary conditions can be imposed either using boundary-adapted basis functions (e.g., Chiavassa & Liandrat 1997) or imposing them afterward, as done for spectral methods using the tau method (Canuto et al. 1987).
We illustrate the wavelet Galerkin method using a Petrov-Galerkin discretization of the 2D incompressible Navier-Stokes equations in vorticity-stream function formulation as proposed by Fr öhlich & Schneider (1996Schneider ( , 1999)).Vorticity-stream function or vorticity-velocity formulations also can be used in the context of wavelet collocation methods (e.g., Vasilyev & Kevlahan 2002).Introducing a classical semi-implicit time discretization of vorticity and stream function equations (here of first order to simplify the presentation) with time step t, we obtain the following system of equations: and The above equations are completed with boundary conditions and a suitable initial condition.At each time step, two elliptic problems (a Helmholtz equation for ω n+1 and a Poisson equation for n+1 ) are solved and a differential operator has to be applied to obtain u n +1 .For spatial discretization, a Petrov-Galerkin scheme is used with orthogonal wavelets as trial functions and operator-adapted wavelets as test functions.The solution is developed into an orthogonal wavelet series (Equation 3), and with the requirement that the residuum vanishes with respect to all test functions, a linear system for the unknown wavelet and scaling function coefficients, ũμ, j i and ū0 i , is solved.The test functions are defined so that the stiffness matrix turns out to be the identity.Therefore, the solution is obtained by simple projection of the nonlinear term.
For the evaluation of the nonlinear term (in which the wavelet coefficients of u n are given), there are two possibilities.The first option is the evaluation in wavelet coefficient space, which requires the computation of connection coefficients,2 as done by Myers et al. (1995), for example.Although many connection coefficients are vanishing or very small, this approach is not efficient in practice (Perrier & Wickerhauser 1999).A second possibility is the evaluation in physical space, similar to pseudospectral methods (Canuto et al. 1987) and hence also called the pseudowavelet technique (Charton & Perrier 1996, Fr öhlich & Schneider 1996, Holmstrom & Walden 1998).Starting from the significant wavelet coefficients of u, the solution is reconstructed on a locally refined grid, and the nonlinear term is then evaluated pointwise.The wavelet coefficients of the nonlinear term can then be computed using adaptive decomposition.
Figure 3 shows an adaptive wavelet computation of decaying homogeneous turbulence in two space dimensions.Comparison of the vorticity with the adaptive grid (locations of active wavelets) demonstrates the automatic grid refinement in regions of strong gradients of vorticity.

Adaptive wavelet collocation methods.
In contrast to adaptive wavelet Galerkin methods that solve problems in wavelet coefficient space, adaptive wavelet collocation methods solve problems in physical space on an adaptive computational grid.This removes two major difficulties associated with wavelet Galerkin methods, namely the straightforward treatment of nonlinearities and general boundary conditions.The evaluation of the nonlinear terms in adaptive wavelet collocation methods is performed in the physical domain, analogous to the pseudowavelet Galerkin methods described above.
Earlier attempts to construct wavelet collocation methods were similar to the wavelet Galerkin approach in that they tried to use operator compression combined with grid adaptation and wavelet interpolation (e.g., Cai & Wang 1996, Prosser & Cant 1998, Vasilyev et al. 1995).In ANRV400-FL42-20 ARI 9 September 2009 21:39 particular, to calculate derivatives, these methods differentiated Approximation 3 and evaluated derivatives of wavelet and scaling functions at collocation points using precomputed values of wavelet derivatives at dyadic grid points.This procedure was proven to be computationally more expensive, especially in multiple dimensions, than finite-difference approximations of the waveletbased interpolant.Many current adaptive wavelet collocation methods take advantage of wavelet interpolation properties and evaluate spatial derivatives using finite-difference approximation on either fixed stencils using a ghost-point approach (e.g., Holmstrom 1999, Vasilyev & Bowman 2000) or variable stencils that use values at the nearest grid points (Rastigejev & Paolucci 2006).Grid adaptation in wavelet collocation methods is done similarly to other wavelet-based methods as discussed above.The only difference is that, at the end of each time step (or iteration), an additional wavelet transform is performed for the analysis of wavelet coefficients and the construction of the adjacent zone.Once Approximation 3 is constructed, the computational grid is obtained based on one-to-one correspondence between wavelets and collocation points: A collocation point is kept if the corresponding wavelet is kept in Approximation 3.
To illustrate the use of wavelet collocation methods in CFD, we discuss the solution of incompressible Navier-Stokes equations in primitive variable form (e.g., Kevlahan & Vasilyev 2005).We note that a similar formulation can be used for the wavelet Galerkin approach (Schneider et al. 2001).To simplify the presentation, we discuss only the first-order time integration of the Navier-Stokes equations.The time-integration scheme is based on a split-step method in which a nonsolenoidal velocity field is first calculated and then is made divergence-free using a pressure projection.In the first step, the nonsolenoidal velocity u(x, t n+1 ) ≡ u n+1 is obtained by solving the linearized equation where L is a linear operator.The Laplacian is discretized implicitly, whereas the advection term is discretized semi-implicitly, which allows the use of a fine grid in the boundary layer, without having to use unreasonably small time steps because of the wavelet-based CFL-like restriction.Once the linear system (Equation 7) is discretized using a wavelet collocation approach (e.g., Prosser 2007, Vasilyev 2003), it can be solved using a variety of methods, e.g., BI-CGSTAB (van der Vorst 1992).In the second step, the velocity u * is corrected by making it divergence-free using the following pressure projection, where P n+1 satisfies The combination of Equations 8 and 9 ensures that the velocity at t n+1 is divergence-free, i.e., ∇ • u n+1 = 0.The Poisson equation (Equation 9) is solved using a wavelet collocation method (e.g., Vasilyev & Kevlahan 2005).Since pressure and velocity are given at the same grid points, one should be careful when choosing the approximation of the Laplace operator.One way to avoid the odd-even decoupling instability associated with nonstaggered grids is to construct the Laplace operator as the inner product of a downwind gradient operator and an upwind divergence operator (Kevlahan & Vasilyev 2005, Schneider et al. 2001).The method was applied for homogeneous turbulent flows in two (Kevlahan et al. 2007) and three (Goldstein et al. 2005) dimensions and for flows around moving and stationary bluff bodies (Kevlahan & Vasilyev 2005).An example of an impulsively started flow through a tightly packed cylinder array at Reynolds number 10 4 is given in Figure 4, where the obstacle was modeled using Brinkman penalization (described in Section 3.6).

Adaptive Multiresolution Methods
Multiresolution-based schemes were originally developed by Harten (1994Harten ( , 1995)), first for hyperbolic conservation laws.Harten's approach was then extended and further developed in different directions by Cohen et al. (2003), Kaibara & Gomes (2000), M üller ( 2003), and Roussel et al. (2003).The main idea of the multiresolution method is to use a hierarchical data representation of the solution.The decay of the wavelet coefficients yields information on the local regularity of the solution.Therewith the truncation error can be estimated, and coarser grids can be used in regions where this error is small and the solution is smooth.An adaptive grid can be introduced by suitable thresholding in which only significant wavelet coefficients are retained.Hence a given discretization on a uniform mesh can be accelerated as the number of costly flux evaluations is significantly reduced, while maintaining the accuracy of the discretization.The memory requirements could also be reduced, for example, by using a dynamic tree data structure.An overview of the different multiresolution methods can be found, e.g., in the books by Cohen (2000) and M üller (2003).Adaptive multiresolution techniques have been mostly employed in the context of compressible flows (for fully adaptive versions with additional memory reduction, the reader is referred to, e.g., M üller 2003).Domingues et al. (2008Domingues et al. ( , 2009) ) developed a fully adaptive multiresolution method using cell averages in combination with finite-volume schemes for the compressible Euler equations, which has recently been extended to the compressible Navier-Stokes equations in two and three space dimensions (Roussel & Schneider 2009).
Finally, we mention that, compared to classical adaptive mesh refinement (AMR) methods (e.g., see Berger & Oliger 1984), which are now widely used for many applications with structured or unstructured grids, multiresolution techniques yield increased data-compression rates in regions where the solution is smooth, thanks to the vanishing polynomial moments of the wavelet functions.Additionally, the grid adaptation error in multiresolution approaches, introduced by thresholding of the details, can be estimated and thusly controlled, which is not the case for AMR methods.

Lagrangian Wavelet Methods
An interesting extension of adaptive wavelet methods is the traveling wavelet method (Basdevant et al. 1990), in which both the position and scale of the wavelets can change continuously in time.Although attractive in principle, the method was proven to be unstable when applied to nonlinear problems, mainly because of wavelet collision.An additional difficulty is associated with the time variability of the scale of the wavelet, resulting in a nonlinear system of ordinary differential equations, even for the case of a linear advection-diffusion problem.Bergdorf & Koumoutsakos (2006) recently overcame this difficulty by enhancing Lagrangian particle methods with multiresolution wavelet-based grid adaptivity.The method extends the adaptive wavelet collocation method (Vasilyev 2003, Vasilyev & Bowman 2000) by combining it with particle methods (Koumoutsakos 2005).This approach leads to a new generation of particle methods with adaptive multiresolution capabilities.A related application for plasma simulations using particle methods for solving the Vlasov equation and wavelet-based grid adaption was proposed by Gutnic et al. (2004).

Space-Time Wavelet Methods
Most of the dynamically adaptive numerical methods, both wavelet-and nonwavelet-based, rely on spatial adaptation, while globally adapting the time step either to ensure stability or to control the time-integration error.This makes these schemes far from optimal for problems that are simultaneously intermittent in both space and time.Different adaptive time-stepping strategies for space-adaptive discretizations of PDEs have been pursued to overcome this difficulty.In the context of adaptive wavelet methods, Bacry et al. (1992) introduced for the first time a scale-dependent time step.They applied this method to linear parabolic equations and the Burgers' equation.For adaptive multiresolution schemes, spatially variable time stepping has been developed to overcome this restriction (Domingues et al. 2008(Domingues et al. , 2009;;M üller & Stiriba 2007) and applied in combination with time-step control (Domingues et al. 2009) for the compressible Euler equations.The underlying idea of these approaches is the following: For explicit time-integration schemes, the time step is imposed by the stability condition on the finest resolution level, whereas in the regions of coarser resolution, the time step can be increased without violating the stability requirement.The missing values at the boundaries between two scales are interpolated in time.
Despite some definite advantages of the local time-stepping methods (especially in the context of adaptive wavelet or multiresolution methods), they are still based on the classical timemarching technique, thus resulting in error accumulation in time, even though the spatial and time-integration errors are controlled at each time step.To address this issue, Alam et al. (2006) proposed a simultaneous space-time adaptive wavelet collocation method in which the problem is solved in a space-time computational domain, which naturally adapts to both space and time resolution to adequately resolve the spatially and temporary intermittent structures of the solution.Besides generating a near-optimal grid for the full space-time solution, this approach also allows the global time-integration error to be controlled.The efficiency and accuracy of the method are demonstrated by applying it to 2D vortex merging (Alam et al. 2006) and 2D homogeneous turbulence (Kevlahan et al. 2007) problems.The space-time method uses substantially fewer (in some instances up to 20 times fewer) space-time grid points and is faster than a comparable adaptive wavelet time-marching method, while achieving similar global accuracy.The main drawback of the simultaneous space-time method is a considerably larger memory requirement, which can be reduced by solving the problem using temporal slices, as described by Alam et al. (2006).

Wavelet-Optimized Adaptive Methods
Both adaptive Galerkin and collocation wavelet methods require the integration of the computational environment, the adaptive multiresolution wavelet transform, wavelet-based interpolation, adaptive grids, and active error control.A computationally efficient implementation of such an integrated approach necessitates the development of a computational environment specifically tailored to adaptive multiresolution wavelet transforms.This is one of the main reasons why many researchers intentionally stayed away from adaptive wavelet methods.However, good waveletcompression properties and the ability to identify localized structures made them an effective tool to be used together with other numerical algorithms.Today, there are a number of waveletoptimized extensions for finite-difference (e.g., Jameson 1998, Leonard et al. 2006), finite-volume (e.g., Liu & Schwarz 2005), and finite-element (e.g., Canuto et al. 1999) methods.In all these extensions, wavelet transforms and analysis of wavelet coefficients are used as a sensor for local mesh refinement.Another wavelet-based acceleration of numerical methods often used in the context of compressible Navier-Stokes solvers is the sparse point representation (SPR) of the fluxes (Chiavassa & Donat 2001), in which fluxes are calculated at a subset of grid points identified by wavelet analysis and interpolated to the rest of the mesh.Different implementations of waveletbased acceleration exist, from simple flux or solution interpolation to more advanced implementations, in which the solution is only advanced at points that belong to the SPR set.The appealing feature of such hybrid approaches is that a wavelet module can be used as an add-on to existing software (e.g., Liu & Schwarz 2005, Sheta et al. 2006).The main disadvantage is that the methodology still requires the use of a fine mesh everywhere in the domain and thus does not reduce memory requirements.Finally, wavelet-based accelerated linear algebra (Beylkin et al. 1991) can be used in the context of any numerical method that requires matrix multiplication or the solution of linear systems.We note, mainly as a precaution, that there are a number of papers in the literature that claim to be wavelet methods, but are based only on single-level approximation using scaling functions (e.g., Qian & Weiss 1993).These methods in no way can be called wavelet methods because they lack the most important aspect, namely multiresolution approximation.Single-level methods, if carefully analyzed, are not that much different from standard finite-difference methods or standard Galerkin methods.

Complex Geometry Treatment
Numerical simulations of flows around solid obstacles of arbitrary complexity are often required for practical engineering applications.There are two general approaches to deal with the complex geometry, namely, body-fitted grid methods (Thompson et al. 1982) and immersed boundary techniques (Mittal & Iaccarino 2005, Peskin 2002).To use adaptive wavelet methods in the context of conventional structured/unstructured body-fitted grid methods, the wavelet transform should be generalized to nontensorial/unstructured meshes that are generated to conform to the complex boundaries.Some progress has been made in this direction, mostly motivated by computer graphics applications (Lounsbery et al. 1997, Reissell 1996, Szczesna 2006, Wang et al. 2008).The use of surface-fitted wavelets for the simulations of flows in complex geometries is mostly limited to SPR of the solution on unstructured meshes (e.g., Valette & Prost 2004).The only use of nontensorial wavelets is the recent extension of the dynamically adaptive wavelet collocation method to secondgeneration spherical wavelets on nested spherical triangular grids (Mehra & Kevlahan 2008).Although the method was only demonstrated for the sphere, it has the flexibility to be extended to other curved manifolds.An alternative approach to simulate flows in complex geometries is to use immersed boundary methods, in which simulations are carried out on nonbody conformal Cartesian meshes, while enforcing complex geometries through penalization.Since Peskin's (1972) immersed boundary method, which was originally introduced to study flow patterns around heart valves, various immersed boundary techniques have been developed.Qian & Weiss (1993) first attempted to combine the single-level wavelet Galerkin method with a capacitance matrix method.Adaptive wavelet Galerkin methods were then extended to solve elliptic and parabolic problems in complex domains using fictitious-domain approaches with boundary conditions imposed through Lagrange multipliers (e.g., Baccou & Liandrat 2006, Diaz 1999, Glowinski et al. 1996, Kunoth 2001).The use of immersed boundary methods in the context of adaptive wavelet methods and Navier-Stokes equations is limited to the Brinkman penalization approach, originally proposed by Arquis & Caltagirone (1984).The main advantage of Brinkman penalization, compared with other immersed boundary methods, is that the error between the solutions of penalized and nonpenalized equations can be estimated rigorously in terms of the penalization parameter (Angot et al. 1999).Furthermore, the solution of the penalized equations converges to the exact solution in the limit as the penalization parameter tends to zero (Angot et al. 1999).Since this volume penalization is simple and cheap to calculate, it is well-suited for flows in complex geometries, including moving or deformable obstacles.The use of adaptive wavelet methods with Brinkman penalization for complex geometry flows was successfully demonstrated for both fixed (Farge & Schneider 2001, Keetels et al. 2007, Schneider & Farge 2002, Vasilyev & Kevlahan 2002, Wirasaet & Paolucci 2005a) and moving (Kevlahan & Vasilyev 2005) obstacles.
Let us briefly outline the main features of the Brinkman penalization method for the incompressible Navier-Stokes equations.Assuming that the problem is defined in a rectangular domain containing obstacles O i , we solve the following set of penalized equations to model the solid boundaries defining O i without explicitly imposing the no-slip condition: with the continuity equation ∇ • u η = 0 and appropriate external boundary conditions.We note that Equation 10 is valid in the entire domain.Here η > 0 is the penalization coefficient, and χ 0 denotes the characteristic (or masque) function that is unity inside the obstacles and zero outside.As η → 0, Carbou & Fabrie (2003) proved theoretically that the error of the solution of the penalized equations (Equation 10) is bounded by u − u η ≤ Cη 1/2 u .Since η is an arbitrary parameter, independent of the spatial or temporal discretization, the boundary conditions can be enforced to any desired accuracy by choosing η appropriately, provided that the penalized equations are properly resolved.This property distinguishes the Brinkman method from other penalization schemes and allows the error to be controlled precisely.
Although it is a flexible and simple method, Brinkman penalization does have a few drawbacks when used with many solvers.First, the large factor 1/η means that the penalty term in Equation 10is stiff and must be solved implicitly, which is easily achieved owing to its diagonal nature.The second drawback is that the penalization term introduces small scales, O((νη) 1/2 ), inside the obstacle in the immediate vicinity of the boundary that needs to be resolved.This often-cited weakness of immersed boundary methods is of concern because of the significant perceived cost associated with the unnecessary resolution introduced outside the boundary region when resolving these small scales in the boundary region.However, this observation is correct only if the penalization method is used together with structured tensorial or zonal meshes, because the overhead in terms of unnecessarily increased resolution outside this region is significant (Mittal & Iaccarino 2005).The use of adaptive wavelet methods reduces this cost to a minimum: The computational overhead is proportional to the volume fraction of the thin skin layer relative to the whole computational volume.The sparseness of the grid inside the obstacle is illustrated in Figure 4b.Moreover, with adaptive wavelet methods, the total simulation error can be maintained to an a priori-prescribed level with optimal grid compression because the automated wavelet-based grid adaptation guarantees the adequate resolution of this skin layer with an a priori-prescribed error controlled by the threshold parameter .The use of immersed boundary methods with adaptive wavelet methods in the context of compressible Navier-Stokes equations is limited to the recent extension of the Brinkman penalization to compressible flows (Boiron et al. 2009, Liu & Vasilyev 2007).To have physically consistent wave reflection from the boundaries, in addition to penalizing momentum and energy equations, the continuity equation for porous media is considered inside obstacles.In the compressible Brinkman penalization approach, the penalized porous region acts as a high-impedance medium, resulting in negligible wave transmissions.

WAVELET ANALYSIS OF TURBULENT FLOWS
A generic feature of turbulent flows is the formation of vortex-like coherent structures whose motions are chaotic, resulting from their mutual nonlinear interactions.To either capture or accurately model vortex dynamics, it is essential to identify and extract coherent structures.Wavelet decompositions of turbulent flows, presented in Section 2, offer an alternative view, as coherent vortices can be defined in an objective way using few wavelet coefficients and without assuming any model for them.Wavelet-based coherent vortex extraction (CVE) allows one to split the flow into two parts: dynamically active coherent vortices, represented by a few wavelet coefficients, and an incoherent background flow, which is noise-like.In addition, the availability of the scale-space information of wavelet coefficients allows one to define statistical tools for the analysis of scale-dependent flow characteristics, such as flow intermittency and anisotropy and local and global energy spectra.These unique features of wavelet-based flow analysis are briefly discussed in this section (for more in-depth discussion, see Abry 1997, Farge 1992, Farge & Schneider 2006, Van den Berg 1999).Farge et al. (1999Farge et al. ( , 2001) ) proposed a wavelet-based method to extract coherent vortices out of both 2D and 3D turbulent flows.This method is based on the vorticity field rather than on the velocity because it preserves Galilean invariance, and simply connected vortex tubes are preserved by nonlinear dynamics due to the Helmholtz theorem.Coherent structures are defined as what remains after denoising.Hence no model is needed for the structures, except for the noise.In the simplest model, the noise is assumed to be additive, Gaussian, and white (i.e., uncorrelated).

Coherent Vortex Extraction
CVE can be viewed as a wavelet denoised vorticity field using Equation 3 with the optimal thresholding parameter opt , which depends on the variance of the incoherent vorticity σ n and on the sample size N.The particular choice of the threshold opt = σ n √ 2 ln N is motivated by denoising theory (Donoho & Johnstone 1994).However, as the variance of the incoherent vorticity is unknown, it has to be estimated from the available total vorticity ω.Azzalini et al. (2005) developed an iterative algorithm based on the method presented by Farge et al. (1999).Typically only one iteration step is performed, which can be justified by the fast convergence of the iterative procedure and by the computational effort being minimized.For notational convenience, the subscript C denotes the coherent component, i.e., ω C = [ω] opt .
The CVE algorithm thus allows one to split the flow into two parts: a coherent flow, corresponding to the coherent vortices, and an incoherent flow, corresponding to the background noise ANRV400-FL42-20 ARI 9 September 2009 21:39 (Farge et al. 1999).This decomposition yields (11) where the subscript I denotes the incoherent component.For orthogonal wavelets, it is easy to show that ω C , ω I = 0 and, hence, the enstrophy is conserved, i.e., Z = Z C + Z I , where Z = 1 2 ω, ω .The corresponding velocity fields can be obtained by inverting the curl operator: The algorithm also works for bi-orthogonal wavelets, with the exception that the enstrophy conservation is not valid anymore.So far CVE has been applied mostly to homogeneous flows in two and three space dimensions.Homogeneous turbulence was decomposed using either orthogonal (Farge et al. 1999(Farge et al. , 2001(Farge et al. , 2003;;Okamoto et al. 2007;Schneider et al. 2006) or bi-orthogonal wavelets (Goldstein & Vasilyev 2004, Roussel et al. 2005).Divergence-free wavelets have also been used (Albukrek et al. 2002, Deriaz & Perrier 2006, Deriaz et al. 2007).Comparisons of CVE with proper orthogonal decomposition-Fourier decomposition were presented by Farge et al. (2003).Jacobitz et al. (2008) extracted coherent vortices in homogeneous shear flows with and without rotation.Temporally developing mixing layers were studied by Schneider et al. (2005), and channel flows were studied by Weller et al. (2006).Lewalle et al. (2000) analyzed experimental mixing-layer data.
To illustrated CVE, we consider its application to 2048 3 DNS data of a 3D homogeneous isotropic turbulent flow (Okamoto et al. 2007).The Reynolds number based on the Taylor microscale corresponds to R λ = 732 (Kaneda et al. 2003).The flow exhibits elongated, distorted, and folded vortex tubes similar to those observed in laboratory (Douady et al. 1991) and numerical experiments (Siggia 1981).Applying CVE yields the coherent and incoherent flow contributions for which isosurfaces of vorticity are shown in Figure 5a.The coherent vorticity is represented by only 2.6% of the wavelet coefficients and retains the same vortex tubes as in the total vorticity, whereas the incoherent vorticity is structureless.The value of the isosurface chosen for visualization is the same for the total and coherent vorticity, but it has been reduced by a factor of 3 for the incoherent vorticity whose fluctuations are much smaller.
The corresponding energy spectra in Figure 5b show that both coherent and incoherent contributions are multiscale.The spectrum of the coherent flow is identical to the one of the total flow all along the inertial range, implying that vortex tubes are responsible for the k −5/3 energy scaling.For the incoherent flow, E(k) is close to k 2 , which corresponds to an equipartition of energy between all wave numbers k.The incoherent velocity is therefore spatially decorrelated, which is consistent with the observation that incoherent vorticity is structureless (Figure 5a).Higherorder statistics of the total and coherent flows also agree well, for both velocity and vorticity, as shown in the corresponding probability density functions in Figure 5c,d.
This example illustrates that CVE works well and yields an efficient representation of coherent vortices with few wavelet coefficients, while preserving the high-order statistics of the flow.In addition, by considering the energy transfer, the nonlinear dynamics is completely preserved all along the inertial range by the coherent contributions (Farge et al. 2003, Okamoto et al. 2007).Okamoto et al. (2007) found that the compression rate (the percentage of wavelet coefficients representing the coherent flow) increases with Reynolds number as N ∝ R 3.9 λ , which proves that the number of modes increases more slowly than for traditional DNS, i.e., N ∝ R 9/2 λ .Goldstein andcolleagues (De Stefano et al. 2005, Goldstein &Vasilyev 2004) studied the effect of the coherent and incoherent contributions on SGS dissipation in the context of LES using a priori and dynamic tests, respectively.Almost all SGS dissipation results from the few SGS modes that make up the coherent part of the SGS field, which indirectly confirms that the coherent component is responsible for the nonlinear flow dynamics.

Scale-Dependent Statistics, Anisotropy, and Relation with Structure Functions
Scale-dependent statistical analysis of turbulent flows can be performed by considering the wavelet coefficients of velocity or vorticity in Equation 1.The coefficients can be interpreted as generalized increments; for example, the classical increments u(t + τ ) − u(t) are recovered by considering the wavelet ψ δ (t) = δ(t + 1) − δ(t).Thus, moments of wavelet coefficients can be directly related to classical structure functions (e.g., Abry 1997, Schneider et al. 2004).Considering ratios of moments, one can define scale-dependent flatness and skewness (Meneveau 1991, Schneider et al. 2004).For orthogonal wavelet decompositions, Parseval's identity holds, which allows the definition of the energy distribution in wavelet coefficient space, depending on scale, position, and direction (Farge 1992, Meneveau 1991).The relation of wavelet spectra to classical Fourier energy spectra has been shown by Farge (1992), Meneveau (1991), andPerrier et al. (1995).Yoshimatsu et al. (2009) studied scale-dependent statistics of fully developed turbulence considering also geometrical interpretations (e.g., the scale-dependent helicity).Bos et al. (2007) developed additional wavelet tools that allow the extraction of directional information.Therefore, one can quantify not only longitudinal and transversal intermittency, but also horizontal or vertical intermittency, which has some interest for anisotropic turbulence, ANRV400-FL42-20 ARI 9 September 2009 21:39 e.g., due to the presence of rotation or stratification.The underlying ideas are straightforward.
Considering moments of one of the three components of the wavelet coefficients ũμ, j i in Equation 1, one can obtain scale and directional information; e.g., second-order moments yield information on the energy distribution as a function of scale and in one of the seven possible directions3 for one of the three velocity components.Ratios of moments allow the definition of scale-dependent flatness and skewness in different directions.Schneider and colleagues (Bos et al. 2007, Yoshimatsu et al. 2009) furthermore demonstrated that the scale-dependent flatness is directly related to the spatial variability of the energy spectrum and thus allows the quantification of the intermittency of turbulent flows.They studied applications for isotropic, rotating, and stratified homogeneous turbulent flows computed by DNS at resolution 512 3 (Bos et al. 2007) and 2048 3 (Yoshimatsu et al. 2009) for isotropic turbulence.

Analysis Based on the Continuous Wavelet Transform
There is a large amount of literature on the analysis of turbulent flows based on the continuous wavelet transform, which is beyond the scope of this review.For more in-depth discussion, the reader is referred to Farge et al. (1990) and Narasimha et al. (2002).Recent work deals with applications of continuous wavelet transforms to vortex-bursting phenomena (Ruppert-Felsot et al. 2009).

HIERARCHY OF TURBULENCE MODELING
Traditional turbulent flow simulations are generally divided into three distinctive classes: (a) fully deterministic simulations, in which all scales of motion are resolved and deterministically computed as in DNS; (b) semideterministic methods, in which some degrees of freedom are computed, while the influence of the others is modeled as in LES; and (c) statistical models, such as RANS methods, in which time-or ensemble-averaged flow fields are computed and the influence of all turbulent fluctuations is modeled.
These three modeling approaches are characterized by different computational costs and a decreasing level of information on the turbulent motions.The most computationally demanding method is DNS, for which the number of degrees of freedom to be computed at each time step is N ∝ Re and N ∝ Re 9/4 for 2D and 3D isotropic turbulent flows, respectively.Therefore, the maximal Reynolds number is limited by the available computational resources (i.e., by memory and CPU performance).The number of degrees of freedom needed by LES, in principle, does not depend on the Reynolds number, as long as one does not resolve the viscous boundary sublayers.Therefore LES allows one to compute much higher-Reynolds number flows than DNS.The least demanding approach in terms of computational cost is RANS, as only the mean flow is computed.
Turbulence is characterized by energetic eddies localized in space and scale, yet the traditional methods discussed above do not take advantage of this localization.In addition to numerical tools for solving the Navier-Stokes equations with self-adaptive discretizations (Section 3), wavelets offer an alternative paradigm to turbulence modeling.The main difference in wavelet-based turbulence modeling is in its compression of the turbulence problem such that a simulation with a subset of the total modes captures the dynamics of the most energetic eddies across the full spectral range of the flow.This section describes a wavelet-based hierarchy of turbulence modeling.The most precise and computationally most demanding approach in this hierarchy is WDNS, which applies the adaptive wavelet-based methods described in Section 3 to solve the Navier-Stokes equations.ANRV400-FL42-20 ARI 9 September 2009 21:39 WDNS (Section 5.1) utilizes the compression properties of the wavelet representation to reduce the computational cost and the memory requirements of DNS, while preserving its accuracy, without performing any modeling.In CVS (Section 5.2), only coherent vortices are resolved.The CVE method described in Section 4 plays a key role in the identification and isolation of coherent vortices.Further reduction of computational cost can be achieved by resolving and tracking only the most energetic part of the coherent eddies, which, similar to LES, necessitates additional modeling of the effect of the less energetic unresolved motions, and results in SCALES (Section 5.3).Finally, we also discuss the application of wavelet-based numerical methods in combination with RANS modeling (Section 5.4).

Wavelet-Based Direct Numerical Simulations
The use of self-adaptive, wavelet-based discretizations of the Navier-Stokes equations described in Section 3 allows one to perform DNS on adaptive grids.The choice of a sufficiently small threshold for wavelet filtering (Equation 3) eliminates the need to do any modeling because the eliminated flow part is not significant for the flow dynamics and the resulting simulations could be interpreted as adaptive DNS (Sagaut et al. 2006).
One of the most important features of the wavelet-based methods, in general, and WDNS, in particular, is the ability to adapt the local resolution to intermittent flow structures and perform simulations with computational complexity considerably lower than commonly accepted.The usual computational estimate of the number of spatial and space-time modes required to simulate decaying 2D and 3D turbulence is N ∝ Re 3/2 and N ∝ Re 3 , respectively.These statistically based estimates neglect intermittency.To determine how N (the number of grid points) and the overall complexity scale with Reynolds number and to verify the claim that WDNS is well-suited for calculating intermittent high-Reynolds number flows, Kevlahan & Vasilyev (2005) conducted a series of 2D simulations of impulsively started flow through a tightly packed cylinder array over a large range of Reynolds numbers, 3 × 10 1 ≤ Re ≤ 10 5 .Figure 4 shows the vorticity and adapted grid at Re = 10 4 .The number of active grid points scales similar to Re 1/2 over five orders of magnitude, whereas the computational complexity scales similar to Re.This represents a significant improvement over the classical estimate of Re 3/2 for 2D turbulence.
Additional studies of 2D decaying turbulence confirmed that the spatial modes scale similar to Re 0.7 , which is slower than the traditional estimate of O(Re) (Kevlahan et al. 2007).A further improvement is obtained by introducing the simultaneous space-time adaptive wavelet method (Alam et al. 2006).Thus the computational complexity of the space-time simulations was N ∼ Re 0.9 for 1.26 × 10 3 ≤ Re ≤ 4.04 × 10 4 over many eddy turnover times (Kevlahan et al. 2007).The relatively high compression confirms the importance of intermittency and encourages further development and application of wavelet-based methods to compute 3D turbulence.

Coherent Vortex Simulation
CVS, introduced by Farge and colleagues (1999;Farge & Schneider 2001) for the simulation of turbulent flows, combines nonlinear approximation with denoising and additionally exploits the properties of wavelets for numerical analysis.As discussed in Sections 3 and 4, wavelets yield attractive discretizations for the Navier-Stokes equations and allow the extraction of coherent vortices in an objective way from turbulent flows.The principles of CVS are (a) to deterministically compute the evolution of the coherent vortices in a wavelet basis, which dynamically adapts to the regions of strong gradients and thus resolves the nonlinear interactions between coherent vortices, and (b) to model the influence of the incoherent components, produced by nonlinear vortex interactions, which are discarded at each time step to model turbulent diffusion.
In the hierarchy of turbulence simulations, CVS falls between DNS and LES.From the filtering viewpoint, what distinguishes CVS from LES is the use of wavelet threshold filters, which depend on the instantaneous flow realization, whereas LES uses linear low-pass filters, which do not adapt to the flow evolution.
The evolution equation of coherent vorticity is obtained from the Navier-Stokes equations by inserting the coherent-incoherent decompositions of vorticity (Equation 11) and velocity (Equation 12) into the equations in vorticity-velocity formulation.Neglecting the influence of the incoherent components ω I and u I , we obtain the following evolution equation for the coherent vorticity ω C : where ν denotes the kinematic viscosity, and ∇ • u C = 0.An orthogonal wavelet threshold filter (Equation 3) with the optimum threshold opt , based on Donoho & Johnstone's (1994) criterion, is applied.The coherent velocity is constructed from the coherent vorticity by applying the Biot-Savart kernel.At the end of each time step, the safety zone described in Section 3 is added and the solution is advanced in time.CVS was successfully applied for 2D and 3D turbulent flows (e.g., Farge & Schneider 2001;Goldstein et al. 2005;Roussel & Schneider 2009;Schneider et al. 1997Schneider et al. , 2006)).
For example, Schneider et al. (2005) compared CVS simulations of a time-developing 3D turbulent mixing layer with DNS, using pseudospectral methods for both computations.For CVS, they filtered the vorticity field each time step and applied the safety zone.This a posteriori test illustrates the potential of CVS without performing fully adaptive computations.
Isosurfaces of the vorticity modulus colored by the spanwise vorticity are plotted in Figure 6.Both visualizations are almost identical; entangled vortex tubes around the two main rolls are present in both DNS and CVS.
The evolution of energy and enstrophy in Figure 6c,d shows that, before the flow is fully turbulent, both quantities are the same for CVS and DNS.After the flow has become fully turbulent, the nonlinear flow dynamics begins to produce incoherent enstrophy.The CVS filtering removes this incoherent enstrophy, explaining why the CVS flow contains less enstrophy than the DNS flow.
Conversely, the time evolution of turbulent kinetic energy E is almost the same for CVS and DNS.Up to the final time, the discrepancy remains less than 0.4% of E. This is a result of the noise-like nature of the incoherent vorticity, which contributes almost nothing to the total energy because it is cancelled out in the Biot-Savart kernel.
The number of wavelet modes retained during CVS (Figure 6e) varies between 8% and 15% N, corresponding to a factor of about two in each Cartesian direction.Approximately 4% N are in the filtered wavelet basis, and the rest correspond to the safety zone, which has been added to track the nonlinear dynamics in space and scale.Figure 6 shows that the number of active wavelet coefficients evolves in time.In particular, when pairings occur, the number of wavelets increases as reflected in the first two peaks in Figure 6e.Roussel & Schneider (2009) have presented fully adaptive CVS computations of a 3D weakly compressible, temporally developing, turbulent mixing layer.An adaptive mulitresolution method based on a second-order finite-volume discretization is used to solve the 3D compressible Navier-Stokes equations in Cartesian geometry (Roussel et al. 2003).CVS were performed by decomposing the conservative flow variables (density, velocity, total energy) into a bi-orthogonal wavelet series by applying the cell-average multiresolution transform.The flow evolution of the coherent part is computed deterministically on a locally refined grid using the adaptive multiresolution method.Figure 7 shows isosurfaces of vorticity together with the corresponding adaptive grid used for the CVS computation.The computational efficiency of CVS in terms of memory and CPU time compression demonstrates that both can be reduced by a factor of three with respect to DNS, whereas energy and enstrophy are reasonably well predicted.

Stochastic Coherent Adaptive Large Eddy Simulations
To further reduce the computational demand of CVS, mainly associated with the requirements of resolving locally coherent modes down to their smallest scale, one has to apply further modeling, which makes the wavelet-based method more computationally feasible for simulations of turbulent flows of engineering interest.In SCALES, introduced by Goldstein & Vasilyev (2004), the formal separation between resolved coherent structures and unresolved eddies is shifted toward the range of more energetic eddies so that the effect of the background flow can no longer be neglected and must be modeled.Introducing an SGS model makes it possible to further reduce the degrees of freedom of the numerical solution, and a higher grid compression with respect to CVS is achieved.In fact, the use of an SGS model makes the SCALES methodology similar to LES.However, in contrast to standard LES based on linear low-pass filters, SCALES exploits a wavelet-based nonlinear thresholding, which depends on the instantaneous flow realization.Furthermore, the distinct difference of the SCALES approach is in the direct coupling of the computational grid and the SGS model: The method has the ability either to compensate for inadequate SGS dissipation provided by the model by increasing the local resolution and hence the level of resolved viscous dissipation or to coarsen the mesh in regions of high SGS dissipation.
The SCALES equations, which describe the space-time evolution of the most energetic coherent eddies in a turbulent flow, can be formally obtained by applying the wavelet thresholding filter (Equation 3) to the incompressible Navier-Stokes equations: where ρ and ν are the constant density and kinematic viscosity, respectively; p is the pressure; and the operator ⊗ denotes the tensor product.As a result of the filtering process, the unresolved quantities, τ = [u ⊗ u] − u ⊗ u , commonly referred to as SGS stresses, are introduced.They represent the effect of unresolved (less energetic) coherent and incoherent eddies on the resolved (energetic) coherent vortices.As in an LES approach, one needs an SGS model to express the unknown stresses in terms of the resolved field to close Equation 15.Similar to LES, the SGS stress has been shown to comprise a minority of coherent modes that dominate the total SGS dissipation and a majority of incoherent modes that, because of their decorrelation with the resolved modes, add little to the total SGS dissipation (De Stefano et al. 2005, Goldstein & Vasilyev 2004).The coherent/incoherent composition of the SGS modes is reflected in the name of the SCALES methodology.Most recent efforts in developing SGS models for SCALES, as for classical LES, have focused on the approximation of the deterministic effect of the coherent SGS motions, whereas the development of stochastic models to simulate the effect of incoherent SGS motions is still an open question.Different SGS models, originally developed in the context of LES, have been successfully extended to the SCALES framework.For instance, there is the dynamic Smagorinsky model, based on the extension of the classical Germano procedure redefined in terms of wavelet thresholding filters (Goldstein et al. 2005).To fully exploit the dynamic adaptivity of the method and to remove the reliance on the presence of homogeneous directions to stabilize the numerical calculations (e.g., Lilly 1992), investigators have recently proposed two classes of local dynamic SGS models for SCALES: a Lagrangian path-line/tube dynamic model (Vasilyev et al. 2008) and local oneequation dynamic models based on the SGS turbulent kinetic energy (De Stefano et al. 2008).The above models make use of the local characteristic filter width, (x, t), which is implicitly defined by the wavelet thresholding procedure (Equation 3) and can be extracted from the analysis of wavelet coefficients during the simulation.It can be interpreted as a measure of the local turbulence resolution.The smaller the value of is, the smaller the turbulence resolution length scale and the greater the fraction of resolved turbulent kinetic energy.In the limit of very small , the WDNS solution is obtained.Such an interpretation of wavelet threshold filtering highlights the similarity between SCALES and classical LES approaches.However, the wavelet filter is different from the LES filters, primarily because it changes in time following the evolution of the solution, which in turn results in an adaptive computational grid that tracks in physical space the areas of dynamically significant flow structures.
Figure 8 summarizes the results of the SCALES simulation of incompressible homogeneous freely decaying turbulence using different SGS models, as well as compares them with CVS, no-model simulations with SCALES threshold, and LES with the global dynamic Smagorinsky model.We note the unique feature of the SCALES approach, namely the coupling of modeled SGS dissipation to grid compression: More grid points are used for models with lower levels of SGS dissipation.In other words, the SCALES approach compensates for inadequate SGS dissipation by increasing the local resolution and, hence, the level of resolved viscous dissipation.This can be seen clearly by comparing the SCALES results using different models with no-model simulations (Figure 8c,d ).Finally, one of the most crucial strengths of the CVS and SCALES approaches is their ability to match the DNS energy and enstrophy density spectra up to the dissipative wave-number range using few degrees of freedom (Figure 8b).Furthermore, the close match for SCALES is achieved using less than 0.4% of the total nonadaptive nodes required for DNS with the same wavelet solver.Comparing these results with those of LES with the global dynamic Smagorinsky model highlights the significance of such a close match.Despite LES using almost four times the number of modes (1.56%), it fails to capture the small-scale features of the spectrum.In addition, with the decrease of SGS dissipation, the spectrum approaches the wavelet-filtered DNS spectrum.Even better approximation of the DNS spectrum is achieved in CVS with the number of modes comparable to LES.

Use of Wavelet-Based Methods with Low-Fidelity Approaches
The use of wavelets with low-fidelity approaches such as RANS is limited to the solution of the modeled equations using adaptive wavelet-based approaches.The computational savings are simply achieved from wavelet compression.The unsteady RANS simulations of the turbulent flow  Stefano et al. 2008, Goldstein et al. 2005, and Vasilyev et al. 2008.around single and tandem cylinders have been performed using both a wavelet-based acceleration approach (Sheta et al. 2006) and the adaptive wavelet collocation method (Liu & Vasilyev 2006).For the latter, the Brinkman penalization method was extended to k − ω unsteady RANS equations with additional penalization terms for turbulent kinetic energy and specific dissipation-rate equations.

CONCLUSIONS AND PERSPECTIVES
Wavelet methods in CFD compose a relatively young area of research.Despite their short, decadelong existence, a substantial number of wavelet techniques have been developed for numerical simulations of compressible and incompressible Euler and Navier-Stokes equations for both inert and reactive flows.What distinguishes wavelet methods from traditional approaches is their ability to unambiguously identify and isolate localized, dynamically dominant flow structures such as shocks, flame fronts, and vortices and to track these structures on adaptive computational meshes.
This article provides a review on the subject of wavelet methods in CFD and discusses different wavelet-based approaches, emphasizing turbulence modeling and simulations using wavelets. 4avelet multiresolution analysis offers a unique framework for modeling and simulating turbulent flows, namely the tight integration of the numerics and physics-based modeling that enables the development of a unified hierarchy of turbulence models of different fidelity.The centerpiece of these models is the energetic coherent structures that capture the dynamics of the flow across the full spectral range.The integration of turbulence modeling with adaptive wavelet methods results in a hierarchical approach, in which all or most energetic parts of coherent eddies are dynamically resolved on self-adaptive computational grids, while modeling the effect of the unresolved incoherent or less energetic modes.Although several open questions still have to be addressed before the methodology can be used for the simulation of high-Reynolds number turbulent flows of engineering interest, it has been demonstrated that adaptive wavelet-based methods for solving Navier-Stokes equations, either in the spirit of adaptive DNS or with wavelet-based turbulence models (such as CVS or SCALES), are more efficient than classical methods in terms of computational complexity with respect to Reynolds number and their ability to capture more flow physics using the same computational resources.The current results found in the literature are promising and show that CPU and memory requirements can be reduced with respect to DNS, while preserving reasonably well the statistical predictability of the computed flows.Nevertheless, for full assessment of wavelet-based turbulence modeling, further detailed benchmark studies are necessary.
Currently we see the largest potential for wavelets in the coupling of adaptive wavelet methods with penalization approaches to compute flows in complex geometries, which may even vary in time or interact with the fluid, as is the case with fluid-structure interactions.The essential properties of wavelet methods, automatic grid generation and error control in combination with turbulence modeling, seem to be most advantageous in these configurations, as dominant flow structures are well localized, even if changing in time, and thus adaptivity pays off.The coupling of multiresolution techniques to accelerate classical solvers is also an interesting direction of great promise.Implementing these multiresolution tools in existing codes is quite straightforward and, similar to AMR techniques, allows the reduction of computational cost, although, in addition to AMR, with an automatic error control.
Adaptive wavelet methods for modeling and simulation of turbulent flows are still in their infancy, and so far mostly flows of academic interest have been studied.Many open questions need to be answered before the methodology can be used as a practical tool for simulations of turbulent flows of industrial relevance.The most pressing issues are the extension of the modeling framework to bounded flows, the definition of coherent structures for highly anisotropic flows, and the efficient implementation of adaptive wavelet solvers on massively parallel computers.In conclusion, we emphasize that wavelet methods in CFD compose a relatively new and evolving area, which involves many challenging issues to be investigated, and we invite researchers to this rapidly advancing field.
Figure 3 Adaptive wavelet computation of decaying 2D turbulence: (a) vorticity with (b) the corresponding adaptive grid.Figure taken from Farge & Schneider 2001.
Figure 4 (a) Vorticity field and (b) corresponding computational grid for the direct numerical simulation of flow around a 2D periodic cylinder array at Re = 10 4 using the adaptive wavelet collocation method.Figure taken from Kevlahan & Vasilyev 2005.
Figure 5 Homogenous isotropic turbulent flow at R λ = 732 and resolution N = 2048 3 .(a) Isosurfaces of total (left panel), coherent (middle panel), and incoherent (right panel) vorticity.The values of the isosurfaces are |ω| = 5σ for the total and coherent vorticity and 5/3σ for the incoherent one with the root mean square σ = (2Z) 1/2 (only subcubes of size N = 256 3 are visualized).(b) Energy spectra of the total, coherent, and incoherent velocity fields.Probability density functions of (c) velocity and (d ) vorticity for the total, coherent, and incoherent components.Figure taken from Okamoto et al. 2007.

Figure 6
Figure 6 Turbulent mixing layer.Isosurfaces of the vorticity modulus colored by the spanwise vorticity at the onset of turbulence: (a) direct numerical simulation (DNS) and (b) coherent vortex simulation (CVS).Time evolution of (c) energy, (d ) enstrophy, and (e) the percentage of retained wavelet coefficients (including the safety zone) used by CVS: DNS ( green line), CVS (red dashed line), and the ratio CVS/DNS (blue dotted line).Figure taken from Schneider et al. 2005.