Wavelet adaptive proper orthogonal decomposition for large-scale flow data

Abstract

The proper orthogonal decomposition (POD) is a powerful classical tool in fluid mechanics used, for instance, for model reduction and extraction of coherent flow features. However, its applicability to high-resolution data, as produced by three-dimensional direct numerical simulations, is limited owing to its computational complexity. Here, we propose a wavelet-based adaptive version of the POD (the wPOD), in order to overcome this limitation. The amount of data to be analyzed is reduced by compressing them using biorthogonal wavelets, yielding a sparse representation while conveniently providing control of the compression error. Numerical analysis shows how the distinct error contributions of wavelet compression and POD truncation can be balanced under certain assumptions, allowing us to efficiently process high-resolution data from three-dimensional simulations of flow problems. Using a synthetic academic test case, we compare our algorithm with the randomized singular value decomposition. Furthermore, we demonstrate the ability of our method analyzing data of a two-dimensional wake flow and a three-dimensional flow generated by a flapping insect computed with direct numerical simulation.

Appendices

Appendix : A. Block-based wavelet adaptation

1.1 A.1 Refinement and coarsening

For the wavelet adaptation scheme, here illustrated for the two-dimensional case, we assume real valued and continuous L²-functions u(x,y), such as the pressure or velocity component of a flow field. The function is sampled on a block-based multiresolution grid as shown in Fig. 1, with maximum tree level ${J}_{\max \limits }$, block size B₁ ×B₂ and a ghost node layer of size g, needed for synchronization. The sampled values on a block ${{{\mathscr{B}}}_{p}^{j}}$ are denoted by:

$$ \begin{array}{@{}rcl@{}} &x^{j}_{p,k_{1}} = {x}_{p} + k_{1} {{\varDelta}} x^{j} &k_{1}=-g,{\dots} ,\text{B}_{1}+g{-1} \end{array} $$

(27)

$$ \begin{array}{@{}rcl@{}} &y^{j}_{p,k_{2}} = y_{p} + k_{2} {{\varDelta}} y^{j} &k_{2}=-g,{\dots} ,\text{B}_{2}+g{-1} \end{array} $$

(28)

$$ \begin{array}{@{}rcl@{}} &\underline{\mathsf{u}}^{j}[p,k_{1},k_{2}] := u\left( x^{j}_{p,k_{1}},y^{j}_{p,k_{2}} \right) \end{array} $$

(29)

For a block refinement j → j + 1 the lattice spacings are halved and dyadic points are added to the block, as shown in Fig. 17.

Fig. 17

Dyadic grid refinement and coarsening of a single block. Refinement: First the block ${{{\mathscr{B}}}_{p}^{j}}$ is refined by midpoint insertion and then split into four new blocks $({{\mathscr{B}}}_{p0}^{j},{{\mathscr{B}}}_{p1}^{j},{{\mathscr{B}}}_{p2}^{j},{{\mathscr{B}}}_{p3}^{j})$. Coarsening: After a low pass filter is applied all midpoints are removed and merged into one block

The values at the refined blocks can be obtained by the refinement relation:

$$ \begin{array}{@{}rcl@{}} \widehat{\underline{\mathsf{u}}}^{j+1} [p,l_{1},l_{2}] = {\sum}_{k_{1}=-g}^{\text{B}_{1}+g{-1}}{\sum}_{k_{2}=-g}^{\text{B}_{2}+g{-1}} h_{l_{1}-2k_{1}} h_{l_{2}-2k_{2}} \underline{\mathsf{u}}^{j} [p,k_{1},k_{2}] \end{array} $$

(30)

where h_k denotes the weights of the one-dimensional interpolation scheme:

$$ \varphi(x) = {\sum}_{k=-N}^{N} h_{k} \varphi(2x-k). $$

(31)

In the following, we will refer to Eq. (30) as the prediction operation $P_{j}^{j+1}:\underline {\mathsf {u}}^{j}\mapsto \widehat {\underline {\mathsf {u}}}^{j+1}$ as it was introduced for point value multiresolution in [28]. Using Eq. (5) one can show that Eq. (4) is equivalent to the continuous refinement relation:

$$ \begin{array}{@{}rcl@{}} \widehat{u}^{(p)}(x,y)& =&{\sum}_{k_{1}=-g}^{\text{B}_{1}+g{-1}}{\sum}_{k_{2}=-g}^{\text{B}_{2}+g{-1}} \underline{\mathsf{u}}^{j}[p,k_{1},k_{2}]\varphi^{j}_{p,k_{1},k_{2}}(x,y) , \end{array} $$

(32)

$$ \begin{array}{@{}rcl@{}} &\text{where}&\quad \varphi^{j}_{p,k_{1},k_{2}}(x,y) = \varphi\left( \frac{x-x_{p,k_{1}}^{j}}{{{\varDelta}} x^{j}}\right)\varphi\left( \frac{y-y_{p,k_{2}}^{j}}{{{\varDelta}} y^{j}}\right) \end{array} $$

(33)

is the two-dimensional tensor product of one-dimensional interpolating scaling functionsφ(x).

When a block is coarsened, the tree level is decimated by one: j + 1 → j and every second grid point is removed. The values at the coarser level are obtained by the coarsening relation:

$$ \begin{array}{@{}rcl@{}} \underline{\mathsf{u}}^{j-1} [p,l_{1},l_{2}] = {\sum}_{k_{1}=-g}^{\text{B}_{1}+g{-1}}{\sum}_{k_{2}=-g}^{\text{B}_{2}+g{-1}} \tilde{h}_{2l_{1}-k_{1}} \tilde{h}_{2l_{2}-k_{2}} \underline{\mathsf{u}}^{j} [p,k_{1},k_{2}] \end{array} $$

(34)

In the notation of Harten [28] coarsening is called decimation and Eq. (34) is denoted by $D_{j}^{j-1}\colon \underline {\mathsf {u}}^{j}\to \underline {\mathsf {u}}^{j-1}$ in the following. After decimation the block will be merged with its neighboring blocks, as shown in Fig. 17. Similar to the continuous refinement relation Eq. (32) there is a continuous counterpart for coarsening: the dual scaling function $\tilde {\varphi }$, which satisfies

$$ \tilde{\varphi}(x) = {\sum}_{k=-N}^{N} \tilde{h}_{k} \tilde{\varphi}(2x-k). $$

(35)

Instead of using Deslauriers-Dubuc (DD) wavelets, as in the framework of Harten [28], we use lifted Deslauriers-Dubuc wavelets, i.e., biorthogonal Cohen-Daubechies-Feauveau wavelets of fourth order (CDF 4,4) with filter coefficients h_k and dual filter coefficients $\tilde {h}_{k}$ listed in Table 1. The lifted DD wavelets allow a better scale separation and can be easily implemented replacing the loose downsampling filter by a low pass filter before coarsening the grid.

Table 1 Filter coefficients h_k and dual filter coefficients $\tilde {h}_{k}$ of the Cohen-Daubechies-Feauveau (CDF 4,4) wavelets applied in the prediction/restriction operation

1.2 A.2 Computing wavelet coefficients and adaptation criterion

With the definition in Eqs. (30) and (34) we have introduced a biorthogonal multiresolution basis $\{\tilde {\varphi }^{j}_{p,k_{1},k_{2}},\varphi ^{j}_{p,k_{1},k_{2}}\}$, which can approximate any continuous function $u\in L^{2}(\mathbb {D})$ arbitrary close. This aspect is the main property of a multiresolution analysis and can be used to relate and compare samples at different resolutions, i.e., different scales.

The difference between two consecutive approximations can be represented by a wavelet with its corresponding coefficients

$$ \begin{array}{@{}rcl@{}} d^{j}[p,k_{1},k_{2}] = \left[u^{j} - P_{j-1}^{j}D_{j}^{j-1} u^{j}\right]_{p,k_{1},k_{2}} = \underline{\mathsf{u}}^{j}[p,k_{1},k_{2}] - \widehat{\underline{\mathsf{u}}}^{j}[p,k_{1},k_{2}] \end{array} $$

(36)

known as wavelet details.

As shown in one space dimension by Unser in [52] for biorthogonal wavelets the difference d^j between two consecutive approximations is bounded for sufficiently smooth L² functions u. The bound depends on the local regularity of u and the order N (here N = 4) of the scaling function:

$$ \|{d^{j}}\| = \|{u - \widehat u}\| \le C_{\varphi,\tilde{\varphi}} ({{\varDelta}} x^{j})^{N} \|{u^{(N)}}\| $$

(37)

where ${{\varDelta }} x^{j}\sim 2^{j}$ is the step size, $C_{\varphi ,\tilde {\varphi }}$ is a constant independent of u. This result can be understood as an interpolation error, since the CDF4,4 filter can be viewed as an interpolation of the averaged data $D_{j}^{j-1}\underline {\mathsf {u}}^{j}$. From Eq. (37) we can thus conclude: since the approximation error of the interpolation scheme depends on the smoothness of the sampled function and the lattice spacing on the block, we can increase the local lattice spacing of the block for blocks where the function is smooth and keep the fine scales only for blocks where u_p is not smooth. This is achieved by coarsening the block, as decreasing the tree level j increases the lattice spacing. For an intended approximation error we therefore define the wavelet threshold𝜖 together with the coarsening indicator $i_{\epsilon }({{{\mathscr{B}}}_{p}^{j}})$:

$$ i_{\epsilon}({{\mathcal{B}}_{p}^{j}}) := \begin{cases} 1 &\text{ if } \underset{k_{2}=0,\dots,\text{B}_{2}{-1}}{\underset{k_{1}=0,\dots,\text{B}_{1}{-1}}{\max}} \mid d^{j}[p,k_{1},k_{2}]\mid<\epsilon\\ 0 & \text{otherwise} \end{cases} $$

(38)

which marks the block for coarsening. For vector-valued quantities $\mathbf{u} =(u_{1},u_{2},\dots ,u_{K})$, a block will be coarsened only if all components indicate coarsening. The pseudocode in Algorithm 1 outlines the wavelet adaptation algorithm for vector-valued quantities. For the sake of completeness we will put this algorithm in relation to the underlying wavelet representation.

1.3 A.3 Wavelet representation in the continuous setting

For completeness of this manuscript we give a detailed description of the underlying wavelet representation in a concise way. For the interested reader we recommend Ref. [50] for a succinct introduction to biorthogonal wavelets and to Refs. [11, 12, 15, 35] for a more detailed description.

Assuming that we have sampled a continuous function $u\in L^{2}(\mathbb {D})$ inside a domain $\mathbb {D}\subset \mathbb {R}^{2}$ on an equidistant grid corresponding to refinement level ${J}_{\max \limits }$, we can block-decompose it in terms of Eq. (29). By choosing $j={J}_{\max \limits }$ in Eq. (32) and summing over all blocks p, we can thus represent u in L² using a basis of dilated and translated scaling functions $\{\varphi _{p,k_{1},k_{2}}^{j}\}$:

$$ u={\sum}_{p\in{\Lambda}^{{J}_{\max}}} u^{(p)}={\sum}_{\lambda\in\overline{\Lambda}^{{J}_{\max}}} c_{\lambda}^{{J}_{\max}} \varphi_{\lambda}^{{J}_{\max}} . $$

(39)

Here we have introduced a multi-index $\lambda =(p,k_{1},k_{2})\in \overline {\Lambda }^{j}:={\Lambda }^{j}\times \{0,\dots ,\text {B}_{1}{-1}\}\times \{0,\dots ,\text {B}_{2}{-1}\}$ for ease of notation. With this notation we can rewrite Eq. (39) in a wavelet series

$$ u ={\sum}_{\lambda\in\overline{\Lambda}^{{J}_{\min}}} c_{\lambda}^{{J}_{\min}} \varphi_{\lambda}^{{J}_{\min}} + {\sum}_{j={J}_{\min}}^{{J}_{\max}-1} {\sum}_{\lambda\in\overline{\Lambda}^{j}}{\sum}_{\mu=1}^{3} d^{j}_{\mu\lambda}\psi_{\mu,\lambda}^{j}, $$

(40)

where the interpolating scaling basis $\{\varphi _{\lambda }^{{J}_{\min \limits }}\}_{\lambda \in \overline {\Lambda }^{{J}_{\min \limits }}}$ approximates u at the coarsest scale ${J}_{\min \limits }$ and the wavelet basis $\{\psi _{\mu ,\lambda }^{j}\}_{\mu =1,2,3,\lambda \in \overline {\Lambda }^{j\ge {J}_{\min \limits }}}$ contains all the additional information necessary to construct u.

In the CDF4,4 setting we have biorthogonal scaling functions:

$$ \begin{array}{@{}rcl@{}} &\varphi_{\lambda}^{j}(x,y)=\varphi\left( \frac{x - x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\varphi\left( \frac{y - y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right), \tilde{\varphi}_{\lambda}^{j}(x,y)=\tilde\varphi\left( \frac{x - x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right) \tilde\varphi\left( \frac{y - y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right)\quad \\ & \quad\text{with}\quad \langle{\varphi^{j}_{\lambda_{1}}},{\tilde{\varphi}^{j}_{\lambda_{2}}}\rangle=\delta_{\lambda_{1},\lambda_{2}} \end{array} $$

and the associated three biorthogonal wavelets (in the d-dimensional case we have 2^d − 1 wavelets see [39, p. 42])

$$ \begin{array}{@{}rcl@{}} &\psi_{1,\lambda}^{j}(x,y)=\psi\left( \frac{x- x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\varphi\left( \frac{y- y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right), \tilde\psi_{1,\lambda}^{j}(x,y)=\tilde\psi\left( \frac{x- x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\tilde\varphi\left( \frac{y- y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right)\\ &\psi_{2,\lambda}^{j}(x,y)=\varphi\left( \frac{x- x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\psi\left( \frac{y- y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right), \tilde\psi_{1,\lambda}^{j}(x,y)=\tilde\varphi\left( \frac{x- x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\tilde\psi\left( \frac{y- y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right)\\ &\psi_{3,\lambda}^{j}(x,y)=\psi\left( \frac{x- x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\psi\left( \frac{y- y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right), \tilde\psi_{1,\lambda}^{j}(x,y)=\tilde\psi\left( \frac{x- x^{j}_{\lambda}}{{{\varDelta}} {x_{p}^{j}}}\right)\tilde\psi\left( \frac{y- y^{j}_{\lambda}}{{{\varDelta}} {y_{p}^{j}}}\right) \\ & \quad\text{with}\quad \langle{\psi^{j_{1}}_{\mu_{1},\lambda_{1}}}, {\tilde{\psi}^{j_{2}}_{\mu_{2},\lambda_{2}}}\rangle= \delta_{\mu_{1},\mu_{2}}\delta_{\lambda_{1}, \lambda_{2}}\delta_{j_{1}, j_{2}} \end{array} $$

for the horizontal (μ = 1), vertical (μ = 2) and diagonal (μ = 3) direction. The wavelet and its dual are defined by the same scaling relations as φ and $\tilde \varphi$:

$$ {\psi}(x) = {\sum}_{k=-N}^{N} {g}_{k} {\psi}(2x-k) \quad \text{and}\quad \tilde{\psi}(x) = {\sum}_{k=-N}^{N} \tilde{g}_{k} \tilde{\psi}(2x-k). $$

(41)

The filter coefficients are given by $g_{k}=(-1)^{1-n}\tilde {h}_{k}$ and $\tilde g_{k}=(-1)^{1-n}{h}_{k}$ with $h_{k},\tilde {h}_{k}$ listed in Table 1. The components of the scaling and wavelet coefficients ($c_{\lambda }^{j}$ and $d_{\lambda }^{j}$) are determined via component-wise projection of Eq. (40) onto the dual basis $\{\tilde \varphi ^{j}_{\lambda }, \tilde \psi ^{j}_{\mu \lambda }\}$:

$$ \begin{array}{@{}rcl@{}} c_{\lambda}^{j} = \langle{u},{\tilde{\varphi}_{\lambda}^{j} }\rangle, \qquad d_{\mu\lambda}^{j} = \langle{u},{\tilde{\psi}_{\mu\lambda}^{j} }\rangle\qquad \text{where} \quad \lambda = (p,k_{1},k_{2}) \end{array} $$

(42)

assuming that $\langle {\varphi ^{j}_{\lambda _{1}}},{\tilde \psi ^{j}_{\lambda _{2}}}\rangle =\langle {\psi ^{j}_{\lambda _{1}}},{\tilde \varphi ^{j}_{\lambda _{2}}}\rangle =0$ are orthogonal. Note that $\langle {a}\rangle {b}={\int \limits }_{\mathbb {D}} a(\mathbf{x})b(\mathbf{x})\mathrm{d} {\mathbf{x}}$ denotes the L²-inner product.

In most wavelet adaptation schemes, one truncates Eq. (40) such that only detail coefficients are kept which carry significant information. According to [45] “this can be expressed as a nonlinear filter”, which acts as a cutoff wavelet coefficients with small magnitude. The cutoff is given by the threshold parameter 𝜖 > 0. However, in contrast to these schemes, our block-based adaptation groups the detail coefficients in blocks, i.e., all details on the block are kept if at least one detail carries important information. This seems to be inefficient at first sight, because unnecessary information is kept, but grouping details in blocks is reasonable, since the block-based adaptation is computationally efficient for MPI distributed architectures. Moreover, groups of significant details are often nearest neighbors, rather than a single significant detail in a block. Therefore, we define the set of blocks

$$ I_{\epsilon}^{j}:=\left \{p \in {\Lambda}^{j}\mid \underset{k_{2}=0,\dots,\text{B}_{2}{-1}}{\underset{k_{1}=0,\dots,\text{B}_{1}{-1}}{\max}}\left| d^{j}[p,k_{1},k_{2}]\right|> \epsilon \right\} $$

(43)

with significant details in the predefined tree level range ${J}_{\min \limits }\le j\le {J}_{\max \limits }$. In the spirit of our previous notation we thus define: $\overline {I}_{\epsilon }^{j}:= I_{\epsilon }^{j}\times \{0,\dots ,\text {B}_{1}{-1}\}\times$ $\{0,\dots ,\text {B}_{2}{-1}\}$ for the set of all significant detail indices. The filtered block-based wavelet field in Eq. (40) can now be written as follows:

$$ u^{\epsilon} ={\sum}_{\lambda\in\overline{\Lambda}^{{J}_{\min}}} c_{\lambda}^{{J}_{\min}} \varphi_{\lambda}^{{J}_{\min}}+ {\sum}_{j={J}_{\min}}^{{J}_{\max}-1} {\sum}_{\lambda\in \overline{I}_{\epsilon}^{j}} {\sum}_{\mu=0}^{3} d^{j}_{\mu,\lambda}\psi_{\mu,\lambda}^{j} $$

(44)

In the following we will denote all fields with an upper index 𝜖, which have been filtered with Algorithm 1 and can be thus expressed as Eq. (44).

For illustration we have computed the vorticity $\omega ^{\epsilon }=\partial _{x} v_{y}^{\epsilon }-\partial _{y} v_{x}^{\epsilon }$ of a thresholded vector field $\mathbf {u}^{\epsilon }=(v_{x}^{\epsilon },v_{y}^{\epsilon },p^{\epsilon })$ in Fig. 3 for various 𝜖 (more details can be found in Section 4.2.1). Here, 𝜖 > 0 and 𝜖 = 0 corresponds to a filtered and unfiltered field, respectively. For increasing 𝜖, less detail coefficients will be above the threshold and therefore the number of blocks decreases.

Taking the difference between the thresholded Eq. (44) and the original field Eq. (10) only details below the threshold are left. Hence the total error can be estimated and yields:

$$ \|{u-u^{\epsilon}}\| \le {\sum}_{j={J}_{\min}}^{{J}_{\max}-1} {\sum}_{\lambda\in {\overline{I}_{\epsilon}^{j}}^{c}} {\sum}_{\mu=0}^{3} \left|d^{j}_{\mu,\lambda}\right|\|{\psi_{\mu,\lambda}^{j}}\| . $$

(45)

Because ${\|\psi _{\mu ,\lambda }^{j}\|}_{\infty }=1$ we finally get ${\| u-u^{\epsilon }\|}_{\infty }\le C\epsilon$ for the total error in the $L^{\infty }$-norm. Similarly one can normalize $\psi _{\mu ,\lambda }^{j}$ in the L²-norm, which corresponds to re-weighting the thresholding criterion $\left |{d^{j}_{\mu ,\lambda }}\right |<\epsilon$ with a level (j) and dimension (d) dependent threshold: $\left |{d^{j}_{\mu ,\lambda }}\right |<\epsilon _{0} 2^{-d(j-{J}_{\max \limits })/2}\epsilon$ [42]. The additional constant 𝜖₀ = 1 for d = 3 and 𝜖₀ = 0.1 for d = 2 tunes the offset of the compression error. It is chosen such that the relative compression error $\mathcal {E}_{\text {wavelet}}$ is close to, but still below 𝜖, i.e., it fulfills Eq. (6) in the L²-norm.

1.4 A.4 L ² Inner products expressed in the wavelet basis

The L² inner product is computed as a weighted sum of two fields u and v. For this we first refine both fields onto a unified grid with identical tree codes Λ^j, as explained in Section 3.2. Then, we are able to compute Eq. (14) as a weighted sum over all blocks:

$$ \begin{array}{@{}rcl@{}} \langle{\mathbf{u}},{\boldsymbol{v}}\rangle &=&{\sum}_{j=1}^{{J}_{\max}}{\sum}_{p\in {\Lambda}^{j}}\langle{\mathbf{u}^{(p)}},{\boldsymbol{v}^{(p)}}\rangle ={\sum}_{j=1}^{{J}_{\max}}{\sum}_{p\in {\Lambda}^{j}}\langle{\boldsymbol{\underline{\mathsf{u}}}^{(p)}},{\underline{\boldsymbol{\mathsf{v}}}^{(p)}}\rangle_{\boldsymbol{I}_{K} \otimes \boldsymbol{W}\otimes\boldsymbol{W}} {{\varDelta}} x_{p}{{\varDelta}} y_{p} \end{array} $$

(46)

$$ \begin{array}{@{}rcl@{}} &&\text{with weights: }(\boldsymbol{W})_{lm} =\langle{\varphi_{l}^{j}},\varphi_{m}^{j}\rangle , \end{array} $$

(47)

Note that this quadrature rule is exact for 𝜖 = 0. We denote by I_K ⊗W ⊗W the Kronecker product between the weight matrix W and the identity matrix $\boldsymbol {I}_{K}\in \mathbb {R}^{K,K}$. The weight matrix is pre-computed by Eq. (47) and its non-vanishing values (W)_ik = w_i−k are shown in Table 2. The listed matrix elements are the discrete values of the autocorrelation function between two compactly supported scaling functions φ, see Fig. 18. Therefore, W is sparse, symmetric and circulant. Since W is also strictly diagonal dominant and all diagonal entries are positive, W and the Kronecker product of such matrices is also positive definite.

Table 2 Values of the autocorrelation function $w_{k}=\int \limits \varphi (x-k)\varphi (x)\d {x}$ of Deslauriers-Dubuc interpolating functions of order N = 2 and N = 4 in the interior of the block. The values for DD4 are rounded to the 4th digit

Fig. 18

Autocorrelation functions of Deslauriers interpolating scaling functions φ of order two (left) and order four (right)

Appendix : B. Derivation of the error estimation given in Eq. (23)

Using Eq. (18) the total error in Eq. (16) becomes

$$ \begin{array}{@{}rcl@{}} \mathcal{E}_{\text{wPOD}}\le \frac{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}} \|{\mathbf{u}_{i} -\mathbf{u}_{i}^{\epsilon}}\|^{2} }{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}} \|{\mathbf{u}_{i}}\|^{2} } + \frac{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}} \|{\mathbf{u}_{i}^{\epsilon} -\tilde{\mathbf{u}}_{i}^{\epsilon}}\|^{2} }{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}}\|{\mathbf{u}_{i}}\|^{2}}. \end{array} $$

(48)

Furthermore we can simplify the first term with Eq. (19) inserting $\|{\mathbf {u}_{i} - \mathbf {u}_{i}^{\epsilon }}\|\le \epsilon \|{\mathbf {u}_{i}}\|$ into the nominator:

$$ \begin{array}{@{}rcl@{}} \frac{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}} \|{\mathbf{u}_{i} - \mathbf{u}_{i}^{\epsilon}}\|^{2} } {{\sum}_{i=1}^{\text{N}_{\mathrm{s}}}\|{\mathbf{u}_{i}}\|^{2}} = \frac{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}}(\epsilon \|{\mathbf{u}_{i}}\|)^{2} }{ {\sum}_{i=1}^{\text{N}_{\mathrm{s}}}\|{\mathbf{u}_{i}}\|^{2} } = \epsilon^{2} . \end{array} $$

The second term in Eq. (48) can be expressed with the help of the eigenvalues of the correlation matrix. We use the identities: ${\sum }_{i=1}^{\text {N}_{\mathrm {s}}}\|{\mathbf {u}_{i}^{\epsilon } - \tilde {\mathbf {u}}_{i}^{\epsilon }}\|^{2}={\sum }_{k=r+1}^{\text {N}_{\mathrm {s}}} \lambda _{k}^{\epsilon }$ for perturbed eigenvalues $\lambda _{k}^{\epsilon }=\lambda _{k}+l_{k}\epsilon$ and ${\sum }_{i=1}^{\text {N}_{\mathrm {s}}}\|{\mathbf {u}_{i}}\|^{2}={\sum }_{k=1}^{\text {N}_{\mathrm {s}}} \lambda _{k}$, yielding

$$ \begin{array}{@{}rcl@{}} \frac{{\sum}_{i=1}^{\text{N}_{\mathrm{s}}}\|{\mathbf{u}_{i}^{\epsilon} - \tilde{\mathbf{u}}_{i}^{\epsilon}}\|^{2}}{{\sum}_{i=1}^{\text{N}_{\mathrm{s}}}\|{\mathbf{u}_{i}}\|^{2}} = \frac{{\sum}_{k=r+1}^{\text{N}_{\mathrm{s}}}(\lambda_{k}+\epsilon l_{k})}{{\sum}_{k=1}^{\text{N}_{\mathrm{s}}}\lambda_{k}} = \mathcal{E}_{\text{POD}}(r,0) + \mathcal{M}_{r} \epsilon . \end{array} $$

(49)

Here ${\mathscr{M}}_{r}={{\sum }_{k=r+1}^{\text {N}_{\mathrm {s}}}l_{k}}/{{\sum }_{k=1}^{\text {N}_{\mathrm {s}}}\lambda _{k}}$ is the perturbation coefficient of the total error. Note that the perturbations l_k are caused by the non-vanishing mixed terms $\langle {\varphi _{\lambda }^{j}}\rangle {\psi _{\mu ,\lambda }^{j}}$ (see also [9]), when computing the correlation matrix from thresholded snapshots $u_{i}^{\epsilon }$. Note that for orthogonal wavelets the first-order perturbations would vanish. For slowly decaying eigenvalues λ_k the perturbation coefficient ${\mathscr{M}}_{r}$ is typically very small, since the sum of perturbations l_k is small compared to the total energy. In this case it is reasonable to neglect the second term in Eq. (49):

$$ \mathcal{E}_{\text{wPOD}}\lesssim \mathcal{E}_{\text{POD}}(r,0) + \epsilon^{2} . $$

(50)

However, in general ${\mathscr{M}}_{r}$ does not vanish and we only have linear convergence in 𝜖:

$$ \mathcal{E}_{\text{wPOD}}\le \mathcal{E}_{\text{POD}}(r,0) + \mathcal{M}_{r}\epsilon+\mathcal{O}(\epsilon^{2}) . $$

(25)

Note that ${\mathscr{M}}_{r}$ does not depend on 𝜖, as all epsilon dependence has been removed. So it is asymptotically a first-order scheme in 𝜖 only. For a certain range we can observe second order, if ${\mathscr{M}}_{r}$ is sufficiently small. Eventually for 𝜖 sufficiently small the first-order term will dominate.

Appendix : C. Technical details and supplementary material

Fig. 19

Example of the light data structure in WABBIT, to be compared with Fig. 2. For each tree lgt_active stores a light-ID list of all active blocks. With the blocks light-ID (lgt_id) all parameters in the forest (tree code, tree ID, tree level, refinement status) can be accessed, from the corresponding row in lgt_block. Note that the order of the light-ID depends on the process rank

Fig. 20

CPU time required for the individual steps of the wPOD algorithm described in Section 3.3. The CPU time is shown for three different cases: 1.) the synthetic test case from Section 4.1 computed with ${J}_{\max \limits }=5$ and 𝜖 = 1.3 × 10^− 4 on a Intel Core i5-7200U cpu (blue), 2.) the flow past a cylinder data Section 4.2.1 computed with 𝜖 = 1.3 × 10^− 4 on Intel Xeon E5645 cpus (orange) and 3.) the bumblebee data of Section 4.2.4 with 𝜖 = 1.3 × 10^− 4 using Intel Xeon Gold 6142 cpus (green). The comparison of the individual cases should be done with care, since they have been computed on different hardware and the data-structure (number of snaphots, blocks, block size, spatial dimension, etc.) is different. Therefore, the computational costs of MPI communication overhead, block administration may vary. However, the overall proportions between the individual steps of the algorithm are comparable among the test cases

Fig. 21

First three sparse modes $\mathbf {\boldsymbol {{\varPsi }}}^{\epsilon }_{k}$ (k = 1 upper left, k = 2 upper right, k = 3 bottom left) with their corresponding amplitudes $a_{k}^{\epsilon }(t_{i})$ computed with 𝜖 = 1.3 × 10^− 2 and one dense (𝜖 = 0) mode $\boldsymbol {{\varPsi }}^{\epsilon }_{3}$ for comparison (bottom, right). Each figure shows the modes “vorticity” (labeled by ∇×v), computed from the two velocity components of the modes, and the pressure component (labeled by p). Note that the first mode represents the base flow, which is non oscillating, whereas the other modes always have an oscillating structure, with a frequency increasing with the mode number. When comparing the dense mode $\boldsymbol {{\varPsi }}^{\epsilon }_{3}$ of the non-adaptive case in the lower right of Fig. 21 with the adapted modes on the lower left, no qualitative differences can be observed, except the local changes in the resolution