Reduced-order modeling (POD)

A design optimization or a parameter sweep solves the same flow problem many times at nearby inputs. DG_windtunnel_model can try a reduced-order model (ROM) before every flow solve. The ROM is a least-squares Petrov-Galerkin (LSPG) solve [CBMF11] on a proper orthogonal decomposition (POD) basis built from the converged solutions computed so far. If the ROM solution is accurate enough, it replaces the full-order (FOM) solve. Otherwise the FOM runs as usual, and its solution is added to the basis.

from dragonfly_sim.core.windtunnel_model import DG_windtunnel_model
from dragonfly_sim.core.reduced_order_model import ReducedOrderModel

rom = ReducedOrderModel(rb_size=20, eta_threshold=1e-10)
model = DG_windtunnel_model(mesh, boundary_dict, ..., reduced_order_model=rom)
model.set_up_sim()
for alpha in alphas:
    u, converged, forces = model.solve_forward(alpha)

The same applies inside a CSDL graph (evaluate), where the snapshot parameters are the shape parameters and the angle of attack. Where the ROM solution is accepted, the derivatives use the adjoint of the reduced equations (see Derivatives below). Elsewhere they use the full-order adjoint.

One evaluation

  1. Once at least min_snapshots (default 2) snapshots exist, build the POD basis \(V\) for the current parameter vector, and project the warm start (the previous accepted solution) onto it.

  2. Solve the LSPG problem in the span of \(V\) (below), and compute the full steady residual of its solution.

  3. Accept the ROM solution if \(\eta = \|R\|/\|u\| <\) eta_threshold.

  4. Otherwise, run the FOM from the warm start exactly as without a ROM. If it converges, its solution becomes a snapshot at the current parameter vector.

With run_fom_every_evaluation=True the FOM also runs when the ROM is accepted, for benchmarking. The output is still the ROM solution. The data store records the ROM wall time, \(\eta\), the ROM coefficients, whether each ROM solve was accepted, the ROM-minus-FOM coefficient errors, and the snapshot singular values (ROM_* and rom_* entries of DataStore).

Snapshots and basis

core/pod.py:SnapshotMatrix keeps the snapshot matrix only as an incrementally updated SVD. It is distributed over MPI with PETSc, and max_rank (default 3 * rb_size) caps the stored rank. A snapshot whose parameter vector repeats an earlier one adds no column. A snapshot that lies in the span of the stored ones, up to a relative residual qr_rank_tol (default \(10^{-10}\), in the inner product below), adds no basis direction; it still counts in the singular values and in the weights of local_weighted. PODBasis builds the basis in one of two ways:

  • basis_mode="global": the leading rb_size left-singular vectors of all snapshots.

  • basis_mode="local_weighted": the leading left-singular vectors of the snapshots weighted by their distance to the current parameter vector. The weights form a partition of unity. The default weight function, "Cubic", is a cubic with compact support, which vanishes at its support radius. By default the support covers the n_nonzero_weights nearest snapshots (default rb_size). With cubic_cutoff=c in \([0, 1)\) the radius is

    \[ r = \max\left(c\,d_\text{min} + (1 - c)\,d_\text{max},\; r_\text{rb}\right), \]

    where \(r_\text{rb}\) is the smallest radius that leaves the rb_size nearest snapshots a nonzero weight: the next distance beyond the rb_size-th smallest, and at least a margin of \(10^{-3}(d_\text{max} - d_\text{min})\) beyond it. The basis therefore always has rb_size columns (as long as the snapshots span that many directions), and a smaller \(c\) lets more snapshots contribute to it. Only snapshots with a basis column count towards rb_size, not repeated parameter vectors. parameter_scales weights the parameters in the distance.

Output-sensitivity weights in an optimization

Shape variables and the angle of attack differ in units and in how strongly the outputs respond to them, so the plain Euclidean distance is dominated by whichever parameter has the largest range. Inside an optimization the adjoint already gives the total derivatives of the objective and constraints, \(\mathrm{d}F/\mathrm{d}p_j\). utils/dv_sensitivity.py reads them without extra solves and ReducedOrderModel.set_distance_sensitivities turns them into one weight per parameter,

\[ s_j = \lVert [\partial F_1/\partial p_j, \dots, \partial F_m/\partial p_j] \rVert_2, \]

normalized to mean 1 and multiplied into the distance (and with parameter_scales). distance_weighting selects 'total_derivative' (default; each output divided by its own magnitude), 'total_derivative_raw' or None. The weights come from the design variables of the CSDL problem, which must be the snapshot parameters (the shape variables, then alpha).

from dragonfly_sim.utils.dv_sensitivity import attach_total_derivative_probe

sim = csdl.experimental.JaxSimulator(recorder, gpu=False)
attach_total_derivative_probe(sim, printer=PETSc.Sys.Print,
                              on_norms=model.reduced_order_model.set_distance_sensitivities)

Without the probe there is no information and the distance is left unweighted, also for basis_mode="local_weighted". The same holds before the first gradient evaluation, and when the weights do not fit the parameter vector. The weights come from the previous design point, so they lag the basis by one evaluation.

Inner product

By default (inner_product="euclidean", state_scaling="freestream"), snapshots and basis use the state-scaled Euclidean inner product

\[ (u, v)_W = \sum_i \frac{u_i v_i}{s_{k(i)}^2}, \qquad W = D^{-2}, \]

with \(s_k\) a typical magnitude of state component \(k\): the free-stream \(\rho\), \(\rho U\) and \(p/(\gamma-1) + \rho U^2/2\), and for RANS the turbulence model’s scale of its variable. Without \(D\), the Spalart-Allmaras variable \(\rho\tilde\nu\) (\(\sim 10^{-4}\) at \(Re = 6\times10^6\)) would carry no weight next to \(\rho E \sim 2\), and a truncated basis would drop it.

inner_product="l2" adds the DG0 mass matrix, \(W = M D^{-2}\), i.e. \(\int_\Omega \sum_k u_k v_k / s_k^2\, d\mathbf{x}\). On NACA0012 angle-of-attack sweeps (Euler at \(M = 0.8\), SA RANS at \(M = 0.7\)) it was 10 to 100 times less accurate than the Euclidean weight. The mass matrix gives the large far-field cells, where the free stream rotates with \(\alpha\), most of the weight, at the expense of the region near the body. state_scaling=None drops \(D\). \(W\) is assembled once, on the undeformed mesh, because the incremental SVD needs a fixed inner product.

LSPG solve

With \(u = V a\) and \(V^T W V = I\), the ROM solves the least-squares Petrov-Galerkin problem \(\min_a \|S R(Va)\|\) by Gauss-Newton:

\[ (S A V)^T (S A V)\, \Delta a = -(S A V)^T S\, r, \qquad a \leftarrow a + \theta\, \Delta a, \]

where \(A\) is the exact Jacobian and \(S = D^{-1}\) (\(D^{-1} M^{-1/2}\) on the current mesh with inner_product="l2") is the residual norm matching the inner product. \(\theta\) is the Newton solver’s positivity-preserving step length, halved until \(\tfrac12\|S r\|^2\) decreases sufficiently (Armijo, line_search=True); if no halving does, the positivity-limited step is taken. The reduced system is small and is solved on every rank.

The stopping rules are those of the original POD benchmark (newton_*, stall_*, max_outer_iterations, g_conv_limit; defaults in brackets):

  • Each pass runs at most newton_max_it [15] iterations. It ends when \(\|g\| = \|(SAV)^T S r\|\) falls below \(\max(\texttt{newton\_atol}, \texttt{newton\_rtol}\,\|g_0\|)\) [1e-12, 1e-10], with \(g_0\) from the first iteration of the first pass (a restart does not reset it); when \(\|\Delta a\| < \max(\texttt{newton\_atol}, \texttt{newton\_rtol} \max(\|a\|, 1))\); or when \(\|g\|\) has not decreased over stall_window [3] iterations (relative stall_rtol [1e-10]).

  • Up to max_outer_iterations [3] passes run, each restarted from the current state. The solve has converged when a pass ends with \(\|g\| < \max(\texttt{newton\_atol}, \texttt{g\_conv\_limit}\,\|g_0\|)\) [1e-10], the same \(g_0\); g_conv_limit = 0 never accepts. Keep g_conv_limit \(\ge\) newton_rtol, or the passes stop before the acceptance target.

\(\|g\|\) is the LSPG stationarity residual. The projected residual \(V^T r\) does not vanish at the LSPG solution, so it is not tested. Acceptance is decided separately, by \(\eta\) against eta_threshold.

Derivatives

The LSPG solution solves the \(k\) projected residual equations

\[ \Phi^T R(V a; x) = 0, \qquad \Phi = S^T S A V, \]

the stationarity condition of \(\min_a \|S R(Va; x)\|\), with \(x\) the shape parameters and the angle of attack. When this solution is accepted, the adjoint of an output \(f\) is the adjoint of these equations, with \(V\), \(S\) and \(\Phi\) frozen at the solution:

\[ (S A V)^T (S A V)\, \lambda = V^T g, \qquad \psi = \Phi \lambda = S^T S A V \lambda, \qquad \frac{df}{dx} = \frac{\partial f}{\partial x} - \psi^T \frac{\partial R}{\partial x}, \]

with \(g = (\partial f/\partial u)^T\). The reduced matrix is the Gauss-Newton matrix of the forward solve. The lifted adjoint \(\psi\) takes the place of the full-order adjoint \(A^{-T} g\), and the rest of the derivative (\(\partial R/\partial \alpha\) and \(\partial R/\partial x_\text{mesh}\), then the mesh warping) is unchanged. \(\psi\) is the Petrov-Galerkin solution of \(A^T \psi = g\) in the span of \(\Phi\), tested with \(V\): \(V^T (A^T \psi - g) = 0\). A Galerkin ROM (\(\Phi = V\)) would give \(V^T A^T V \lambda = V^T g\) and \(\psi = V \lambda\) instead. DG_windtunnel_model.apply_inverse_jacobian computes \(\psi = S Q y\) from a distributed thin QR \(S A V = Q R\) and \(R^T y = V^T g\), so it never forms the normal matrix, whose condition number is squared. This costs one Jacobian assembly and \(k\) products with it, and no Krylov solve. CSDL passes one right-hand side per output, and every output after the first reuses the factorization. DataStore.adjoint_solves records which adjoint each derivative used.

This is the derivative of the reduced equations with the basis held fixed. Two kinds of terms are left out:

  • The exact linearization of the stationarity condition also has second-derivative terms of \(R\) contracted with \(z = S^T S R(u_r)\), the weighted residual at the ROM solution: \(V^T \big(\sum_i z_i \nabla^2_{uu} R_i\big) V\) in the reduced matrix, and \(\partial_x \big(z^T A(u; x) V\lambda\big)\) in the derivative. With inner_product="l2", \(S\) also depends on the mesh. These terms are linear in the residual of the accepted solution, so they are of relative size \(O(\eta)\).

  • For basis_mode="local_weighted" the basis depends on \(x\) through the snapshot weights, and that dependence is not differentiated. For "global" the basis is fixed.

CSDL re-runs the model before every derivative evaluation, so each design is evaluated at least twice. A re-evaluation of the design that the full-order model solved last skips the ROM. The full-order solve restarts from its own converged solution, so it takes no Newton steps, and the derivative is the full-order adjoint. By then that solution is a snapshot, and the ROM would reproduce it and be accepted. A design the ROM solved keeps the ROM and its reduced adjoint.

What to expect

The ROM pays off when the solutions over the parameter range lie close to a low-dimensional linear subspace (a small Kolmogorov \(n\)-width). Shocks that move with the parameters are the hard case, the more so at p = 0, where they are smeared over a few cells. Angle-of-attack sweeps over \(1.4^\circ\) to \(2.8^\circ\) on the NACA0012 meshes in meshes/, with the default settings (rb_size=20), 10 evaluations and the full-order solve run every time for comparison, gave:

Case

\(\eta\) of the ROM, 2 → 9 snapshots

\(c_l\) error

wall time, ROM / full order

Euler, \(M = 0.8\), 240,752 dofs

\(2.4\times10^{-5}\) → \(2.0\times10^{-6}\)

\(2\times10^{-2}\) → \(10^{-4}\) to \(10^{-3}\)

1–5 s / 26–84 s

SA-neg RANS, \(M = 0.7\), \(Re = 6\times10^6\), 47,400 dofs

\(1.9\times10^{-5}\) → \(1.1\times10^{-7}\)

\(10^{-2}\) → \(10^{-5}\) to \(10^{-3}\)

1–4 s / 17–40 s

All ROM solves converged, but none met the default eta_threshold of \(10^{-10}\). On these cases every evaluation fell back to the full-order solve, and the ROM added a few seconds per evaluation. Raise eta_threshold only where the coefficient errors at that \(\eta\) are acceptable; run_fom_every_evaluation=True measures them.

Flow models

The ROM uses only the host interface of CompressibleEulerModel, which CompressibleRANSModel forwards:

u_vec

the state, whose layout the basis rows follow

problem.F, problem.J

residual and exact Jacobian; for RANS the condensed-gradient problem, which updates the reconstructed gradient

steady_residual()

switches the pseudo-time (PTC) and physical time terms off for the reduced solve

positivity_step_length(x, dx)

the Newton step limiter’s step length

state_scales()

the magnitudes \(s_k\) (a new flow model with extra variables supplies them through the extra_scales plug-in)

There is therefore no model-specific ROM code. The reduced solve is plain damped Gauss-Newton, also for the RANS model, whose full-order steady solve uses pseudo-transient continuation. This is a sound choice because the ROM starts from the projected previous solution.