3.1. Characterizing the quasi-periodic dynamics

To probe the quasi-periodic dynamic behaviour, we work within a spatially discrete formulation of the governing equations, which we write generically as

(3.1)\begin{equation} \boldsymbol{\mathsf{M}}\dot{\boldsymbol{q}}=\boldsymbol{f}(\boldsymbol{q}). \end{equation}

In (3.1), $\boldsymbol {q}$ conveys the flow state sampled at a subset of points $n_s$ within the flow domain $\varOmega$; e.g. in an incompressible setting, $\boldsymbol {q}$ is a concatenation of the spatially discrete flow state variables. (For an incompressible flow, these would be the various velocity components and pressure at the collection of discrete spatial points used to represent $\varOmega$.) The vector $\boldsymbol {f}$ is the dynamical evolution operator. The matrix $\boldsymbol{\mathsf{M}}$ allows for non-unity weightings on the time derivatives, including a singular $\boldsymbol{\mathsf{M}}$ as is needed to represent conservation of mass in the incompressible Navier–Stokes equations (resulting in a differential–algebraic system of Hessenberg index 2).

We assume that the dynamics (3.1) contains underlying near-periodic dynamic behaviour, and that the user has a means of identifying distinct cycles of this behaviour. That is, IPhaB POD requires that an operation $\boldsymbol {\mathcal {F}}$ be known such that $\boldsymbol {q}_l=\boldsymbol {\mathcal {F}}(\boldsymbol {q})$ provides a meaningful way to probe the intrinsic near-periodic content (without the additional dynamics driven by this near-periodic content). That is, $\boldsymbol {q}_l$ exhibits the targeted near-periodic dynamics that allows for distinct cycles, and the phase along each cycle, to be identified. We refer to the dynamics of $\boldsymbol {q}_l$ (i.e. the dynamics of $\boldsymbol {q}$ under $\boldsymbol {\mathcal {F}}$) as restricted dynamics. While the requirement of an appropriate $\boldsymbol {\mathcal {F}}$ may appear somewhat limiting, in many cases such knowledge is available. Examples of flows with identifiable near-periodic content include wakes, separation bubble ‘breathing’, fluid–structure interaction problems, and rotating systems. We will next characterize mathematically the assumed near-periodic behaviour, and give examples of this mathematical representation for the two representative cases described in § 2.

We assume that the restricted dynamics $\boldsymbol {q}_l=\boldsymbol {\mathcal {F}}(\boldsymbol {q})$ yields $n_c$ cycles. We further assume that the restricted trajectory associated with the $k$th cycle, $\{\boldsymbol {q}_l(t): t\in [t_{0_k}, t_{f_k} ]\}$, $k=1, 2, \ldots, n_c$, has a (near-)period $T_k:= t_{f_k}- t_{0_k}$. Associated with each $k$th cycle is also a minimum and maximum at, say, $t_{min_k}$ and $t_{max_k}$, given by derivative conditions of $\boldsymbol {q}_l$ over that cycle. The peak-to-peak amplitude for each cycle can thus be defined as $A_k := \|\boldsymbol {q}_l(t_{max_k})-\boldsymbol {q}_l(t_{min_k})\|$. One can define a nominal period duration and peak-to-peak amplitude from the arithmetic mean of the quantities $T_k$ and $A_k$ over all cycles, e.g. $\bar {A}= (1/n_{c})\sum _{k=1}^{n_c} A_k$ (with an analogous definition for $\bar {T}$).

Because the dynamics is assumed to be imperfectly periodic, in general $\boldsymbol {q}_l(t_{f_k})\ne \boldsymbol {q}_l(t_{0_k})$. In fact, there may be cycles for which the period duration is very different from the nominal period $\bar {T}$, and the restricted trajectory $\{\boldsymbol {q}_l(t): t\in [t_{0_k}, t_{f_k} ]\}$ yields large values of $\|\boldsymbol {q}_l(t_{f_k})- \boldsymbol {q}_l(t_{0_k})\|$. For the purposes of this paper, we focus on restricted trajectories that exhibit near-periodic orbits, and omit extreme events from consideration. That is, we perform IPhaB POD only on cycles that are near-periodic in the sense that their restricted trajectories in phase space are closed, and their associated time durations are equal to a representative period duration (to within some tolerance). The near-periodic trajectories considered within IPhaB POD are then computed as those satisfying

(3.2)\begin{equation} \gamma_{j} = \{\boldsymbol{q}_l(t): t \in [t_{0_{j}},t_{f_{j}}],|\bar{T}-T_j|/\bar{T} := \delta_j<\delta^*,\|\boldsymbol{q}_{l}(t_{f_{j}}) - \boldsymbol{q}_{l}(t_{0_{j}})\|/\bar{A} := \epsilon_j <\epsilon^*\}, \end{equation}

where $\delta ^*$ and $\epsilon ^*$ are threshold values for departure from the nominal period duration and from trajectory closedness. We denote the number of near-periodic orbits satisfying (3.2) as $n_p$, the set of indices associated with these orbits as $N$ (with $\text {dim}(N)=n_p\le n_c$), and $\varGamma =\{\gamma _j\}$, $j\in N$.

It will be useful in presenting IPhaB POD to define the ‘most’ representative periodic cycle of the restricted dynamics, which we denote as the trajectory $\gamma ^p := \{\boldsymbol {q}_l(t): t\in [t_{0_k}, t_{f_k} ]\}$, $k\in N$, that satisfies $\min _{k} \sqrt {\delta _k^2 + \epsilon _k^2}$. (Note that by construction, $\gamma ^p$ is one of the admitted near-periodic trajectories in (3.2).) We denote the period and discrepancy from trajectory closedness of the nominally periodic $\gamma ^p$ as $T^p$ and $\epsilon ^p$, respectively.

We now discuss how the trajectories $\gamma _j$, $j\in N$, and $\gamma ^p$ are obtained for the two problems considered in this paper. In addition to illustrating how these trajectories are determined, the examples are also intended to motivate as modestly reasonable the assumption that a means to probe the near-periodic content within the dynamics, given by a user-provided mapping $\boldsymbol {\mathcal {F}}$, is available.

For the aerodynamic flat plate problem, vortex shedding is the (only) content within the dynamics, and is perfectly periodic to within numerical error (cf. figures 1(a–d) and the surrounding discussion). The lift dynamics can be used to characterize this vortex-shedding behaviour. That is, while the lift is not equivalent to the large-scale vortex shedding, its dynamics can be used to identify distinct phases of this vortex shedding (i.e. periodic) behaviour. Thus, for this problem, an $\boldsymbol {\mathcal {F}}$ can be defined that takes the flow state, extracts the flow stresses via the gradient of the pressure field and the viscous Laplacian of the flow velocity, and evaluates these on the surface of the aerodynamic body. (Or alternatively, many immersed boundary methods provide the surface stresses directly as a state variable, and under this formulation the surface stress state variable may be spatially integrated.)

The periodic trajectories $\varGamma$ can be extracted from this definition of $\boldsymbol {\mathcal {F}}$ as follows (as illustrated in figure 2). The lift signal depicted in figure 2(a), plotted sufficiently far after the initial transient that the limit cycle has been reached, demonstrates the periodic behaviour of the system. The variation in trajectory closedness of the lift cycle and discrepancy from the nominal period duration are shown in figure 2(b). For this problem, the cycle-to-cycle differences are small and due to discretization errors, so thresholds $\epsilon ^* = 5 \times 10^{-5}$ and $\delta ^* = 2 \times 10^{-3}$ are chosen to include all cycles. For illustration purposes, and to indicate a phase trajectory undergone by the lift as a representation periodic dynamics, figure 2(c) shows the various lift cycles superimposed on one another as a phase portrait of ${\rm d}C_l/{\rm d}t$ versus $C_l$, though only a single curve is visible because of how nearly periodic the system is. The cycle depicted in red in figures 2(a,c) is the one computed as $\gamma ^p$ based on the criterion described above, though again the periodicity of the system dynamics renders the specific choice of $\gamma ^p$ non-essential for this problem. Each distinct near-periodic trajectory within $\varGamma$, defined through the above procedure, indicates a distinct vortex-shedding cycle. Moreover, the trajectory exhibited by each $\gamma _j\in \varGamma$ experiences each of the different phases of the vortex-shedding process within a single cycle.

Figure 2. (a) The coefficient of lift $C_l$ for the representative cycle of vortex shedding from the flat plate is shown in red, and the remaining ensembles are plotted in green. (b) Nearness of cycle start and end positions ($\epsilon$) and difference of cycle period to nominal period ($\delta$) for all cycles. Threshold values used to determine near-periodic cycles are shown as dashed lines (for this case, all cycles are considered nearly periodic). (c) Cycles plotted using a phase portrait with the coefficient of lift plotted against its time derivative. All cycles lie on top of each other because all are nearly periodic.

Note that at higher Reynolds numbers, the vortex shedding of flow past a flat plate would lose perfect periodicity due to turbulent fluctuations, though it would retain a strong signature at the characteristic vortex-shedding frequency that is robust to Reynolds number when scaled appropriately (Roshko Reference Roshko1954). For this higher Reynolds number case, one may wish to modify the $\boldsymbol {\mathcal {F}}$ operation to not only compute the lift but also remove the high-frequency components associated with the turbulent fluctuations. This low-frequency signal could be used to establish the (near-)periodic portion of the vortex-shedding behaviour that distinguishes the various cycles.

For the shock problem, the (near-)periodic behaviour is the undulating motion of the shock, with aforementioned finer-scale features such as the supersonic jet exhibiting dynamics about this near-periodic behaviour (cf. figures 1e–h). For this problem, our choice of $\boldsymbol {\mathcal {F}}$ leverages an image filter: an edge-detection algorithm is used to define the upstream (downstream) region of the shock, labelling this region with value zero (unity). An unaltered schlieren image is shown in figure 3(a), with the filtered image using the edge-detection algorithm for that same image shown in figure 3(b). The filtered image is then integrated in space to yield a one-dimensional signal that varies in time, shown in figure 3(c). (The time instance shown in figures 3(a,b) is indicated by the red marker in figure 3(c).) This spatially integrated signal encodes the near-periodic portion of the dynamics, as its value changes as the shock moves within the domain. Similar to the flat plate problem, to indicate phase trajectories undergone by the integrated signal from figure 3(c), figure 3(e) depicts the integrated signal against its time derivative.

Figure 3. (a) Schlieren image of shock oscillation over cone-cylinder body. (b) Filtered version of (a) to identify phase of near-periodic dynamics. (c) Spatially integrated value of the filtered image shown in (b) plotted against time, with the time instance of (a) marked in red. (d) Nearness of cycle start and end positions ($\epsilon$) and difference of cycle period to nominal period ($\delta$) for all cycles. Threshold values used to determine near-periodic cycles are shown as dashed lines. (e) Cycles plotted using a phase portrait with the value of plot (c) plotted against its time derivative. Cycles that are considered near-periodic using the thresholds of (d) are shown in green, while the representative cycle is shown in red.

Visible in both figures 3(c,e) is an apparent lack of periodicity that nonetheless contains common cyclic behaviour (with significant cycle-to-cycle variation). This cycle-to-cycle variation is demonstrated quantitatively in figure 3(d) through the difference between each cycle's start and end values ($\epsilon$) and the difference between each cycle's period and a reference period ($\delta$). For this problem, the thresholds of these values that are used to define near-periodic orbits were chosen to be $\epsilon ^* = 0.1$ and $\delta ^*= 0.05$.

These thresholds were chosen as a balance between (1) removing extreme events with significant differences from the nominal trajectory based on closedness and period length, and (2) over-trimming cycles with some differences due to cycle-to-cycle variations that are important to understand. The characteristic periodic cycle $\gamma ^p$ is highlighted in red in figure 3(e), while all cycles are shown in black, and all cycles that are the near-periodic cycles belonging to $\varGamma$ are shown in green. The computed trajectories $\varGamma$ contain information about the shock motion for cycles where the undulations remain near a representative periodic orbit. As such, each restricted trajectory $\gamma _j\in \varGamma$ represents a distinct cycle of the shock's oscillation process, and the trajectory exhibited by each of these cycles reflects the various shock locations explored through an oscillation cycle.

Note that the data that we utilize for the shock problem come from schlieren imaging, so that the filtering operation defined above is performed not on the flow state $\boldsymbol {q}$ directly, but on integrated observations of refractive index changes induced by density gradients in the flow. We do not focus here on this distinction; one could envision an analogous $\boldsymbol {\mathcal {F}}$ operation that acts on the full compressible flow state, identifying discontinuities in any of the flow quantities, and using this location to label upstream and downstream regions and integrate the result over the spatial domain as indicated above.

We highlight that while in both examples the operation $\boldsymbol {\mathcal {F}}$ is chosen to map from the flow state to a real scalar value, this choice is not unique (nor is the choice of mapping, even when restricting the image of the operation to the collection of real numbers). Indeed, mappings onto higher-dimensional outputs might facilitate a more robust analysis of the near-periodic behaviour and the activity driven by this behaviour, described below. As the main focus of this paper is on utilizing a relevant representation of the intrinsic near-periodic dynamics on which to perform POD, we leave the pursuit of an optimal way to extract the underlying near-periodic content as a point for future research.

3.2. A POD formulated for near-periodic driving behaviour

This subsection describes how to use the collection of near-periodic trajectories $\varGamma$, associated with the restricted dynamics, to define an expectation operator and inner product that compare distinct cycles along matched phases of the underlying near-periodic content. The result is an energy-ranked modal decomposition, optimal in the usual sense of POD of the best rank-$k$ approximation, for any $k\in \mathbb {Z}^+$, for a data set organized in the proposed phase-matched way.

The expectation operator is defined as an ensemble average over the collection of full-flow-state trajectories that, under the map $\boldsymbol {\mathcal {F}}$, belong to the near-periodic trajectories $\varGamma$. That is, the sample space being drawn from is $j\in N$, and we may, for clarity, explicitly write the full flow state as a random variable or a function of this sample space via

(3.3)\begin{equation} \tilde{\boldsymbol{q}}(t, j) = \boldsymbol{q}(t), \quad \text{for } t\in[t_{0_j}, t_{f_j}],\ j\in N, \end{equation}

where the dependence on the sample space is stated explicitly.

Intuitively, we seek an expectation operator defined approximately as an ensemble average over these realizations $\langle \tilde {\boldsymbol {q}} \rangle (t) = ({1}/{n_p}) \sum _{j\in N} \tilde {\boldsymbol {q}}(t, j)$. (Recall that $n_p=\text {dim}(N)$.) One issue with this form is that it is not possible to evaluate this expectation operation at a given time instance, since each $\tilde {\boldsymbol {q}}(\boldsymbol {\cdot}, j)$ is defined only over the time interval $t\in [t_{0_j}, t_{f_j}]$. More substantively, even barring that challenge, comparing across the same time within a cycle does not generally correspond to the same phase of the near-periodic dynamics, because of the deviations from periodicity of the nominally periodic behaviour.

To address both of these issues, we define a prototypical periodic trajectory, and relate the time interval of each of the cases with near-periodic cycles $j\in N$ to the nearest time instance of the reference prototypical cycle. This characteristic trajectory is chosen as $\gamma ^p$, defined in § 3.1, for which the associated $\sqrt {\delta ^2+\epsilon ^2}$ is minimized (i.e. minimizing the difference from the nominal period length and the departure from trajectory closedness).

To enable a sensible averaging operation expressed in terms of a unified time interval, we require a map from each time instance $t\in [t_{0_j}, t_{f_j}]$ to the time instance along the nominal trajectory, $t^p\in [t_{0}^p, t_{f}^p]$, for which the near-periodic dynamics are most similar. Denoting $\tau _j = [t_{0_j}, t_{f_j}]$ and $\tau ^p = [t_{0}^p, t_{f}^p]$, we define this map $\mathcal {T}_j^p: \tau _{j} \rightarrow \tau ^p$ using

(3.4)\begin{equation} t^p = \min_{t' \in [t^p_0, t_f^p]} \| \boldsymbol{q}_l^p(t') - \boldsymbol{q}_{l_{j}}(t) \|, \quad t\in\tau_j. \end{equation}

There is a related map $\mathcal {T}_p^j: \tau ^{p} \rightarrow \tau _j$, where $\mathcal {T}_j^p = (\mathcal {T}_p^j)^{-1}$ provided that $\mathcal {T}_p^j$ is injective. Intuitively, this requirement is modest: the trajectories are selected to represent nearly the same phenomenon, so supposing a one-to-one mapping from one time interval to the other is natural.

Using (3.4), we arrive at an expression for the expectation operator, given by

(3.5)\begin{equation} \langle \tilde{\boldsymbol{q}} \rangle(t) = \frac{1}{n_p} \sum_{j\in N} \tilde{\boldsymbol{q}}(\mathcal{T}_j^p(t), j), \quad t \in \tau^p. \end{equation}

This expression mirrors a phase average, with the phase of the near-periodic dynamics identified using (3.4). Usually, when phase averages are performed, phase information is known a priori through some external forcing such as the rotation of an aerofoil or a rotor. In this case, the phase is identified from the analysis of the instantaneous snapshots: it is this ‘intrinsic phase’ identification and definition that inspired the name IPhaB POD. We will refer to the outcome of (3.5) as the IPhaB mean, to distinguish it from the time mean and from more common phase averages that are based on an externally driven phase variation. (The terminology also fits the name of the method, as the expression (3.5) is formally the mean under the expectation operator used in our proposed IPhaB method.) This definition of an expectation operator is also strongly connected to the conditional expectation defined in Schmidt & Schmid (Reference Schmidt and Schmid2019). Indeed, the expression (3.5) entails a conditional expectation operator. Unlike in Schmidt & Schmid (Reference Schmidt and Schmid2019), however, the expectation here is conditioned not on a specific event, but rather on a phase of an intrinsic near-periodic behaviour. This difference has required new formulations to the definition of the realizations developed in (3.2) as well as a time mapping that facilitates a clear comparison of the dynamics of each realization relative to the underlying near-periodic behaviour. This time mapping directly affects the definition of the expectation operator in (3.5), and will also manifest itself in the definition of the inner product (described next).

With the expectation operator defined, POD can be performed once a suitable inner product is chosen. Similar to the expectation operator, we emphasize that an inner product that simply integrates along the trajectory of each realization, $\tilde {\boldsymbol {q}}(\boldsymbol {\cdot}, j)$ ($\,j\in N$), would not account for the imperfect periodicity of the underlying dynamics, and would yield a modal representation that erroneously uses ‘higher modes’ to convey this underlying cyclical behaviour. We therefore choose to perform this trajectory integration relative to a nominal periodic trajectory. Each instance of the trajectory encoded in $\tilde {\boldsymbol {q}}(\boldsymbol {\cdot}, j)$ ($\,j\in N$) will be cast in terms of the nearest instance on the nominal trajectory. This approach facilitates a meaningful comparison across the various ensembles: all trajectory integrations are aligned with a shared (near-)orbit reflective of the periodic dynamics that each of the ensembles is close to, as defined by (3.2). Denoting the manifold on which the dynamics $\boldsymbol {q}$ evolves as $\boldsymbol {\mathcal {M}}$, we define the inner product as

(3.6)\begin{equation} (\boldsymbol{r},\boldsymbol{m}) = \int_{t_0^p}^{t_{f}^p} \boldsymbol{m}^*(\mathcal{T}_j^p(t'))\,\boldsymbol{\mathsf{W}}\, \boldsymbol{r}(\mathcal{T}_j^p(t')) \,{\rm d}t', \quad \forall \boldsymbol{r}, \boldsymbol{m}\in \boldsymbol{\mathcal{M}}, \end{equation}

where $(\boldsymbol{\cdot})^*$ denotes complex conjugation, and $\boldsymbol{\mathsf{W}}$ is an appropriately chosen positive-definite weighting matrix that incorporates quadrature weights to approximate spatial integration of the space-discrete state variables, and allows for the induced norm to be a meaningful representation of physical energy. (Because of the space–time formulation employed here, most applications of this method would likely involve real-valued functions. We leave the complex conjugate in the inner product definition for generality.) The exact construction of $\boldsymbol{\mathsf{W}}$ varies from problem to problem, but is generally dependent on the field value chosen. The state variable used for the flat plate problem is vorticity, so the weighting chosen is a diagonal matrix with ${\mathsf{W}}_{jj}=\Delta x_j\,\Delta y_j$, the space between the $j$th spatial point and its adjacent one. That is, the resulting norm represents the enstrophy (per unit of density) of the system. The ‘state variable’ used in the shock problem is the pixel intensity provided by each schlieren image. The image intensity is driven by density gradients in the flow, and the same weight ${\mathsf{W}}_{jj}=\Delta x_j\,\Delta y_j$ is chosen. (Note that the third spatial dimension is neglected for both the two-dimensional flat plate problem and for the three-dimensional shock problem, since for the latter, only planewise data are available.) This choice avoids weighting certain regions of space more than others. While the resulting norm has a non-obvious direct physical interpretation because the pixel intensity is not a flow state variable, it is proportional to the square density gradient field.

The eigenvalue problem that yields the IPhaB POD modes then takes the standard form, obtained by searching for the mode $\boldsymbol {\phi }$ that maximizes the expectation of the projection of the flow state onto that mode. The result is

(3.7) \begin{equation} (\mathcal{R}\boldsymbol{\phi})(t) = \langle (\boldsymbol{\phi}, \tilde{\boldsymbol{q}}')\,\tilde{\boldsymbol{q}}' \rangle (t) = \int_{t_0^p}^{t_f^p} \boldsymbol{\mathsf{R}}(t, t')\,\boldsymbol{\mathsf{W}}\,\boldsymbol{\phi}(t') \,{\rm d}t' = \lambda\,\boldsymbol{\phi}(t), \quad t \in \tau^p, \end{equation}

where, as is common practice, the eigenvalue problem is defined with respect to the fluctuation state variable $\tilde {\boldsymbol {q}}' = \tilde {\boldsymbol {q}} - \langle \tilde {\boldsymbol {q}} \rangle$. Note also that $\boldsymbol{\mathsf{R}}(t, t') = \langle \tilde {\boldsymbol {q}}'\tilde {\boldsymbol {q}}'^* \rangle (t, t')$ is the autocorrelation function.

Because the proposed formulation utilizes the standard POD formulation for a specific choice of inner product and expectation operator, prior results of operator self-adjointness, optimality and mode coherence apply. The operator $\mathcal {R}$ is self-adjoint, and the eigenvalue problem (3.7) yields a collection of modes $\{\boldsymbol {\phi }_1, \boldsymbol {\phi }_2, \ldots \}$ that are orthogonal with respect to the inner product (3.6). The associated eigenvalues are positive and real-valued, and are taken here to be ordered such that $\lambda _1 \ge \lambda _2 \ge \cdots$ (see e.g. Holmes et al. (Reference Holmes, Lumley, Berkooz and Rowley2012) for a review of these results). An equation for, say, the $n+1$ most energetic modes can be obtained by maximizing the expectation of the projection of the flow state onto a subspace orthogonal to the span of the first $n$ modes. Optimality of low rank representation and space–time coherence of the modes also hold. For a discussion of the latter point about coherence, and the breakdown of this property when inner products that do not incorporate time are used, see Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018). Note that in our case, coherence is with respect to time scaled by the nominal near-periodic trajectory.

3.3. Computing IPhaB POD modes with sampled data

To approximate the time-continuous formulation above with discrete samples of data, we construct a data matrix where the $j$th column contains the flow state associated with the quasi-periodic trajectory $\gamma _j$, $j\in N$. We recall that to facilitate a meaningful comparison across the other trajectories, the eigenvalue problem (3.7) is cast in time according to the nominal periodic trajectory $\gamma ^p$ rather than the $j$th trajectory. Denoting by $n_n$ the number of time samples of this nominal periodic trajectory, and recalling that $n_p=\text {dim}(N)$ is the number of near-periodic cycles contained within the set $\varGamma$, the data matrix to be formed is then

(3.8)\begin{equation} \boldsymbol{\mathsf{X}}=\begin{bmatrix} \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_1}^p(t_1), n_1) & \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_2}^p(t_1), n_2) & \ldots & \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_p}^p(t_1), n_p)\\ \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_1}^p(t_2), n_1) & \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_2}^p(t_2), n_2) & \ldots & \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_p}^p(t_2), n_p)\\ \vdots & \vdots & \ddots & \vdots \\ \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_1}^p(t_{n_n}), n_1) & \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_2}^p(t_{n_n}), n_2) & \ldots & \tilde{\boldsymbol{q}}'(\mathcal{T}_{n_p}^p(t_{n_n}), n_p) \end{bmatrix} . \end{equation}

Note that each block entry in the matrix $\boldsymbol{\mathsf{X}}$ is defined with respect to the fluctuating quantity. For example, for cycle $n_j\in N$ and time instance $t_k\in \tau _p$, the rows are built as $\tilde {\boldsymbol {q}}'(\mathcal {T}_{n_j}^p(t_{k}), n_j) = \tilde {\boldsymbol {q}}(\mathcal {T}_{n_j}^p(t_{k}), n_j) - \langle \tilde {\boldsymbol {q}} \rangle (\mathcal {T}_{n_j}^p(t_{k}))$.

Using the definition (3.8) of the data matrix, the eigenvalue problem that approximates (3.7) is

(3.9)\begin{equation} \boldsymbol{\mathsf{X}}\boldsymbol{\mathsf{X}}^*\tilde{\boldsymbol{\mathsf{W}}}\boldsymbol{\phi}_d = \boldsymbol{\varLambda} \boldsymbol{\phi}_d, \end{equation}

where $\tilde {\boldsymbol{\mathsf{W}}}$ is a weighting matrix that incorporates the matrix $\boldsymbol{\mathsf{W}}$ from (3.7) as well as quadrature weights to approximate the time integral in the inner product definition. This weighting matrix is block diagonal, and constructed as a concatenation of $({1}/{n_p})\boldsymbol{\mathsf{W}}$, repeated along block diagonal entries $n_p$ times. This scaling accounts for the ensemble average (and the spatial integral via $\boldsymbol{\mathsf{W}}$). (For reference, for results shown in later sections for the commonly used space-only POD, which forms the data matrix to have each snapshot of the full data set in a separate column, the weighting matrix is $\tilde {\boldsymbol{\mathsf{W}}} = ({1}/{n_t})\boldsymbol{\mathsf{W}}$, where $n_t$ is the number of snapshots used for the space-only decomposition. This number need not equal $n_n\times n_p$, since generally space-only POD does not omit snapshots from different cycles, and instead uses the full data set.) The eigenvector $\boldsymbol {\phi }_d$ discretely approximates the time-continuous eigenvector $\boldsymbol {\phi }$. With the state vector $\boldsymbol {q}\in \mathbb {R}^{n_s}$, where $n_s$ is the number of points used to represent the spatial domain $\varOmega$, $\boldsymbol {\phi }_d\in \mathbb {R}^{n_s\times n_n}$. The $k$th collection ($k=1, \ldots, n_n$) of $n_s$ entries in $\boldsymbol {\phi }_d$ approximates $\boldsymbol {\phi }(t_k)$. In practice, the eigenvalue problem (3.9) is typically computationally unwieldy, and the method of snapshots (Sirovich Reference Sirovich1987) may be used to extract the desired modes from the smaller eigenvalue problem associated with $\boldsymbol{\mathsf{X}}^*\tilde {\boldsymbol{\mathsf{W}}}\boldsymbol{\mathsf{X}}$.

We now detail the process used to compute $\mathcal {T}_{j}^p(t_k)$, where $j=n_1, \ldots, n_p$, $k=1, \ldots, n_n$. First, we define from the data the nominal near-periodic trajectory $\gamma ^p=\{\tilde {\boldsymbol {q}}(t_{0_l}, l), \tilde {\boldsymbol {q}}(t_{1_l}, l), \ldots, \tilde {\boldsymbol {q}}(t_{f_l}, l)\}$ for the index $l\in N$ that yields the minimal value of $\sqrt {\epsilon ^2+\delta ^2}$. (Note that this trajectory $\gamma ^p$ is the time-sampled variant of the trajectory defined in § 3.2; we avoid using a distinct variable for notational ease.) That is, the time instances indicated in the data matrix (3.8) are defined as $\{t_0, t_1, \ldots, t_{n_n}\}:=\{t_{0_l}, t_{1_l}, \ldots, t_{f_l}\}$. We then build the data matrix where each column contains a different realization of $j\in N$, $j\ne l$.

For each $j$th column, we identify $\mathcal {T}_{j}^p(t_k)$ at each time instance $t_k$ along the nominal near-periodic cycle as follows. First, we identify for the $j$th cycle the time $t_q$ that yields a restricted flow state nearest to the restricted flow state associated with the nominal periodic cycle:

(3.10)\begin{equation} t_q = \operatorname{arg\,max}_{t'\in\tau^j} \frac{\boldsymbol{\mathcal{F}}(\tilde{\boldsymbol{q}})(t_{k}, l)^*\, \boldsymbol{\mathcal{F}}(\tilde{\boldsymbol{q}})(t', j)}{\|\boldsymbol{\mathcal{F}}(\tilde{\boldsymbol{q}})(t_{k}, l)\|^2}. \end{equation}

For some physical problems, including the cone-cylinder shock problem studied in this work, the approach (3.10) can struggle when the near-periodic trajectory of the restricted dynamics contains multiple points in phase space that are close together or overlapping within a single cycle. In the case of the shock problem, this issue occurs when the shock passes the same location, with a similar overall shape, during the upstream and downstream phases of the shock motion, respectively. In this work, this challenge was managed by further restricting the set of times $t'\in \tau ^j$ considered for the optimizer $t_q$ in (3.10) to be within a given range from the prior identified cycle time $t_{q-1}$. This additional restriction prevents erroneous matching of phases at opposite sides of the dynamic cycle. Future work will aim to resolve this challenge more robustly through a dynamical systems perspective.

Because of finite sampling in time, this nearest flow state in the $j$th cycle to that associated with the desired phase of the near-periodic dynamics will not in general be desirably close to $\boldsymbol {\mathcal {F}}(\tilde {\boldsymbol {q}})(t_{k}, l)$. As such, a locally quadratic representation in time, about $t_q$, is used to find the $t\in [t_{q-1}, t_{q+1}]$ with the maximal projection onto $\boldsymbol {\mathcal {F}}(\tilde {\boldsymbol {q}})(t_{k},l)$. Once this time instance $\mathcal {T}_{j}^p(t_k)$ is determined, quadratic interpolation is used to populate $\tilde {\boldsymbol {q}}'(\mathcal {T}_{j}^p(t_k),j)$ from the values of the fluctuation flow state $\tilde {\boldsymbol {q}}'(t_{q-1})$, $\tilde {\boldsymbol {q}}'(t_q)$ and $\tilde {\boldsymbol {q}}'(t_{q+1})$. This process is used to construct $\boldsymbol{\mathsf{X}}$ in (3.8), for $j=n_1, \ldots, n_p$, $k=1, \ldots, n_n$.

This proposed approach yields, at each cycle (column in the data matrix $\boldsymbol{\mathsf{X}}$), snapshots that are matched according to the phase of an underlying near-periodic behaviour driving the dynamics. This approach is useful for problems where the driving behaviour is near-periodic, as it enables a meaningful comparison across cycles whose strong similarities (a common near-periodic behaviour) are obscured by detailed differences. The proposed approach also has utility in cases where the dynamics is perfectly periodic. Even in this setting, finite sampling in time will generally yield a collection of snapshots within a cycle that are at different instances along the phase of the periodic trajectory, which if compared directly do not enable a meaningful assessment of how the dynamics evolve relative to this periodic behaviour.