AeroSim Docs
Mathematical foundations of PODFS¶
The PODFS method rests on two pieces of mathematics that are standard in data analysis but are not part of the everyday toolkit of every CFD practitioner: the proper orthogonal decomposition, which is an eigenvalue problem in disguise, and the Fourier series. This page builds both from the ground up, assuming only a working familiarity with vectors and integrals. It is a primer, not a derivation of the method itself; the method and its use in Nassu are on the main PODFS page.
A snapshot is a vector¶
The starting point is to stop thinking of a velocity field on the inlet plane as a picture and start thinking of it as a single vector. Suppose the inlet plane is discretised into \(P\) points and at each point we store the three fluctuating velocity components. Stacking all of those numbers into one long column gives a vector with \(3P\) entries,
one such vector per recorded instant \(t_{m}\). A precursor recording of \(M\) snapshots is therefore \(M\) points in a space of dimension \(3P\) - a very large dimension, but the points do not fill it. Turbulence is organised: the snapshots cluster near a low-dimensional surface inside that space, and the whole purpose of POD is to find the few directions along which they actually vary and discard the rest. Everything below is a way of making “direction” and “few” precise.
Inner products, norms and orthogonality¶
To talk about directions we need the geometry of that space, and the geometry comes from the inner product. For two snapshot vectors \(\mathbf{a}\) and \(\mathbf{b}\) it is the familiar sum of products of corresponding entries, which for fields is written as an integral over the plane,
From it follow the two notions the method needs. The norm \(\lVert\mathbf{a}\rVert = \sqrt{\left\langle\mathbf{a},\mathbf{a}\right\rangle}\) measures the size of a field; for a velocity fluctuation \(\lVert\mathbf{u}'\rVert^{2}\) is twice its kinetic energy integrated over the plane, which is why POD energy talk is exact and not a metaphor. Orthogonality is the statement \(\left\langle\mathbf{a},\mathbf{b}\right\rangle = 0\): two fields are orthogonal when they have no overlap, no component of one lying along the other. A set of fields that are mutually orthogonal and each of unit norm is called orthonormal.
Why orthonormal bases are the convenient ones
A basis is just a set of building-block fields from which any snapshot can be assembled as a weighted sum. If the basis is orthonormal, finding the weights is trivial: the weight of building block \(\boldsymbol{\phi}_{n}\) in a field \(\mathbf{u}'\) is simply the inner product \(\left\langle\mathbf{u}',\boldsymbol{\phi}_{n}\right\rangle\) - a projection, exactly like reading off the \(x\)-component of an arrow by projecting it onto the \(x\)-axis. There is no system of equations to solve, the weights are independent of one another, and the energy of the field is the plain sum of the squared weights (Parseval’s relation). A non-orthogonal basis loses every one of these properties, which is the deeper reason the main page warns that modes from two different decompositions must never be blended: they are not mutually orthogonal, so their weights stop being independent and the energy bookkeeping breaks.
Eigenvalues and eigenvectors¶
The second idea is the eigenvalue problem. A square matrix \(\mathbf{C}\) acting on a vector usually rotates and stretches it. For special vectors, however, the action is pure stretching with no rotation: the output points the same way as the input. Such a vector is an eigenvector \(\boldsymbol{\psi}\), and the stretch factor is its eigenvalue \(\lambda\),
The eigenvectors are the natural axes of the matrix - the directions it treats simply - and the eigenvalues say how strongly it acts along each. When \(\mathbf{C}\) is symmetric (equal to its own transpose) two facts hold that POD relies on completely: the eigenvalues are all real and non-negative, and eigenvectors belonging to distinct eigenvalues are automatically orthogonal. So a symmetric matrix hands us, for free, an orthogonal set of directions ranked by a real, non-negative number. That is precisely the orthonormal-basis-plus-energy-ranking that the previous section said we wanted.
The covariance picture
The matrix POD diagonalises is a correlation (covariance) matrix, and for a covariance matrix the eigen-decomposition has a vivid meaning. Picture the cloud of snapshot points; it forms an elongated, tilted ellipsoid. The eigenvectors are the principal axes of that ellipsoid - the direction of greatest spread, then the greatest spread orthogonal to it, and so on - and each eigenvalue is the variance (the energy) along its axis. Ordering the eigenvalues from large to small orders the axes from “where the data varies most” to “where it barely varies at all”. Keeping only the first few axes is keeping the few directions along which the turbulence actually lives. This is identical to principal component analysis in statistics; POD is that same procedure applied to a velocity field.
Proper orthogonal decomposition¶
POD is the marriage of the two ideas. It seeks the orthonormal basis \(\left\{\boldsymbol{\phi}_{n}\right\}\) - the modes - that captures the most energy in the fewest terms, and the answer turns out to be the eigenvectors of the two-point correlation of the fluctuation field. Each snapshot is then written as a weighted sum of modes,
where the spatial modes \(\boldsymbol{\phi}_{n}\) are fixed shapes and the temporal coefficients \(a_{n}\left(t\right)\) are the projections that say how strongly each shape is present at each instant. The eigenvalue \(\lambda_{n}\) is the energy carried by mode \(n\), so ordering the modes by eigenvalue orders them by energy, and truncating after \(N\) modes keeps the \(N\) most energetic shapes - the optimal spatial compression in the sense that no other \(N\)-term linear basis captures more energy Lumley[1], Berkooz et al.[2].
In practice the correlation acts on the \(3P\)-dimensional snapshot space, which is enormous. The method of snapshots Sirovich[3] sidesteps this: because there are only \(M\) snapshots, the energetic directions all lie in the span of those \(M\) fields, and the modes can be obtained from the much smaller \(M\times M\) matrix of snapshot-to-snapshot inner products,
This matrix is symmetric and positive semi-definite by construction, so it has the real non-negative eigenvalues and orthogonal eigenvectors promised above. Its eigenvectors give the temporal coefficients directly, and the spatial modes are recovered by combining the snapshots with those eigenvector weights. Diagonalising an \(M\times M\) matrix (a few thousand snapshots at most) instead of a \(3P\times3P\) one (millions of grid entries) is what makes POD tractable for a CFD inlet plane.
Fourier series¶
POD compresses space. The Fourier series compresses time. Its premise is that any reasonably smooth signal that repeats with period \(T\) can be written as a sum of sines and cosines whose frequencies are integer multiples of the base frequency \(1/T\),
where the integer \(k\) counts how many full oscillations fit in one period and the complex coefficient \(c_{k}\) sets the amplitude and phase of that harmonic. Two properties make this the natural representation for a periodic turbulent coefficient. First, the harmonics are themselves orthogonal over the period - the same orthogonality idea as for the spatial modes, now applied to functions of time - so each coefficient \(c_{k}\) is found by an independent projection of the signal onto its harmonic and the harmonics do not interfere. Second, a smooth signal’s coefficients shrink rapidly with \(k\), so keeping only the lowest \(2K+1\) harmonics reproduces the signal to high accuracy while collapsing a long sampled time series to a short list of numbers. That truncation is the temporal compression.
What the Fourier step buys, concretely
Storing each temporal coefficient as a Fourier series rather than as raw samples is what makes the replayed inflow continuous (the series is a closed-form function, evaluable at any instant, not a table that must be interpolated between stored steps), periodic (it loops seamlessly with period \(T\) by construction), and time-step-independent (the downstream solver picks its own \(\Delta t\) and simply evaluates the series there) Treleaven and others[4]. For a lattice-Boltzmann solver, whose time step is fixed by the lattice speed of sound and generally differs from the precursor’s, this removes the need to resample stored data onto the lattice clock every iteration - the single strongest reason PODFS suits an LBM inlet.
Putting the two together¶
PODFS is the composition of the two compressions, applied to the fluctuation only. POD turns a long sequence of full-plane snapshots into a handful of spatial modes \(\boldsymbol{\phi}_{n}\left(\mathbf{x}\right)\) and their temporal coefficients \(a_{n}\left(t\right)\); the Fourier step turns each of those coefficients into a few harmonics \(c_{n,k}\). Reassembling them gives the inlet field as the mean plus a Fourier-weighted sum of modes,
which is exactly equation (6) on the main page. Read with the geometry of this primer in mind: the modes are an orthonormal, energy-ranked basis for the inlet plane, each mode’s strength in time is a short Fourier series, and the whole field is rebuilt by projecting back out of that doubly-compressed representation. Two orthogonal decompositions, one in space and one in time, are all the method is.