Physics-Informed Neural Networks

This page documents the neural-network half of the project: how the surrogate is designed, how it is trained on solver output, and the engineering iterations that took it from a toy model to a parametric design-space surrogate. Every design decision is justified by a diagnostic measurement, in the same spirit as the C++ solver documentation.

1. Component Reference

2. Network Architecture & Design

The surrogate is a fully-connected multilayer perceptron. Three architectures share the same Fourier-feature backbone, growing from a single operating point to a continuous spatio-temporal design space:

# Steady-state PINN (single flow condition)
Input:  (x, y)                       normalized to [-1, 1]
Hidden: 2 -> 64  (x8 layers, tanh) -> 3
Output: (u, v, p)

# Parametric PINN (continuous design space)
Input:  (x, y, Re_n)                 Re_n = (Re - 100) / 900
Fourier: frozen random sinusoidal projection (m=128, sigma=5.0)
       lifts (x, y) -> 512-dim frequency space
Hidden: 513 -> 256 (x8 layers, tanh) -> 3   (462K params)
Output: (u, v, p)

# Time-parametric PINN (spatio-temporal surrogate)
Input:  (x, y, Re_n, t_n)            t_n = frame / (n_frames - 1)
Fourier: frozen random sinusoidal projection (m=128, sigma=5.0)
       lifts (x, y) -> 512-dim frequency space
Hidden: 514 -> 256 (x8 layers, tanh) -> 3   (593,155 params)
Output: (u, v, p)

# Time-parametric PINN -- multi-scale Fourier upgrade (current)
Fourier: 3 frozen bands (sigma=1, 5, 20, m=128 each)
       lifts (x, y) -> 768-dim frequency space
Hidden: 770 -> 256 (x8 layers, tanh) -> 3   (856,067 params)
Loss:   NS + pressure-Poisson + vorticity + data + BC + IC
Train:  1k pretrain + 8k Adam + 1k L-BFGS, adaptive sensors

The training objective is a hybrid loss that fuses data and physics:

L_total = w_data * MSE(pred, LBM)
        + w_pde  * NS_residual(pred)
        + w_bc   * BC_loss(pred)

The PDE residual is the steady incompressible Navier-Stokes equations, differentiated with torch.autograd. This is what separates a PINN from a plain data-driven surrogate: the network is penalized whenever its field violates conservation of mass and momentum, so it generalizes into regions where solver data is sparse. Boundary-condition losses are case-specific (cavity: no-slip walls + moving lid; cylinder: inflow, outlet, walls, and surface). Typical weights are w_pde = 1.0, w_data = 10.0, w_bc = 5.0.

3. Training Pipeline

A training run proceeds as follows:

1. load frame*.json          NumPy fields (u, v, p) on the downsampled grid
2. sample sensors            importance-sampled points near the geometry
3. sample collocation        PDE residual points across the interior
4. normalize                 (x, y) -> [-1, 1], Re -> Re_n, frame -> t_n
5. optimize                  Adam (lr=1e-3) + CosineAnnealingLR
6. evaluate                  full-grid inference -> 3-panel comparison

For the cavity case, training is multi-Re: each epoch samples collocation points at Re = 100, 400, and 1000, so a single model learns the whole slice of the design space rather than one operating point. The time-parametric variant extends this to the full transient: collocation points are sampled across all 51 frames of the LBM evolution, the PDE loss enforces the unsteady Navier-Stokes equations (with material time derivatives), and an initial-condition loss pins the rest state at t_n = 0. A data-only pretraining phase runs first to avoid the constant-collapse failure mode, where the network learns to output the zero-mean field that simultaneously minimizes both the data loss and the PDE residual. Sensors are adaptively resampled by residual every 2,000 epochs to track the moving vortex core. On the cylinder case, 40% of collocation points are placed within three cylinder radii via importance sampling to resolve the near-field gradients. All training runs on the Apple M5 using the PyTorch MPS (Metal Performance Shaders) backend; the cylinder steady-state run completes in roughly 9 minutes for 15,000 epochs, while the cavity time-parametric run takes about 223 minutes (1,000-epoch pretrain + 8,000 Adam + 1,000 L-BFGS steps).

# Install dependencies (Python >= 3.10)
pip3 install -r requirements.txt

# Smoke-test the data loader (NumPy only)
python3 data/loader.py

# Train the cylinder steady-state surrogate
python3 cases/cylinder/train.py

# Train the cavity parametric surrogate (multi-Re)
python3 cases/cavity/train_steady.py

4. Engineering Iterations: Design & Fixes

The cavity parametric surrogate was built through three iterations. Each step states the hypothesis, the change made, the metric that moved, and the diagnosis that motivated the next step.

Iteration 1 -- 64-wide MLP, no Fourier features

Hypothesis: a small MLP can interpolate the cavity field across Re = 100 and 400 directly from normalized coordinates.

Change: parametric PINN with input (x, y, Re_n), 8 hidden layers of 64 units (116K parameters), raw coordinates (no Fourier window).

Result: u L2 of 73.5% at Re = 100 and 52.6% at Re = 400; v L2 45.6% / 38.8%; p L2 107% / 108%.

Diagnosis: the network is simultaneously capacity-limited and spectrally biased. It captures the broad recirculation but cannot represent the lid-driven shear layer or the pressure variation, and the error is worst at the moving-lid boundary.

Iteration 2 -- 256-wide MLP, no Fourier features

Hypothesis: the L2 errors are a capacity problem; widening the network will resolve the gradients.

Change: 8 hidden layers of 256 units (462K parameters), still raw coordinates.

Result: u L2 improved to 56.1% / 41.8%; v L2 34.7% / 33.0%; p L2 12.6% / 12.6%.

Diagnosis: more capacity helps the pressure field but the velocity field is still wrong in character. The predicted peak velocity is only 39% of the true value (0.038 vs 0.097), and the predicted pressure span is just 1.7% of the true span, with error concentrated at the lid boundary layer. A fine-tune with L-BFGS improved the loss by only 0.2% -- so this is a representation problem, not an optimization problem. The tanh MLP simply cannot fit the high-frequency lid shear.

Iteration 3 -- Fourier feature window

Hypothesis: spectral bias is the root cause; feeding the network high-frequency coordinate features breaks it.

Change: a frozen random sinusoidal projection (m = 128 basis functions, sigma = 5.0) lifts (x, y) into a 512-dimensional frequency space before the MLP. The parametric input becomes 513-dimensional (512 Fourier + Re_n); the MLP stays at 256 x 8.

Result: measured u L2 of 23.7% at Re = 100 and 24.4% at Re = 400; v L2 29.3% / 30.0%; p L2 12.5% / 12.6%. The u_max ratio (pred/true) dropped from 3.50 (iteration 2) to 1.24 (Re=100) / 1.10 (Re=400), confirming the spectral bias is removed.

Diagnosis: the Fourier window gives the optimizer high-frequency basis functions to combine, so the lid shear layer and the vortex core are no longer averaged away. This is the architecture carried forward as the project baseline.

Iteration 4 -- Time-parametric (spatio-temporal) extension

Hypothesis: a single network can learn the full transient, not just a steady frame, delivering an ML-powered animation that exceeds the LBM section in interactivity while keeping the same spatial fidelity.

Change: extend the input with a normalized time t_n = frame / (n_frames - 1), giving a 514-dimensional input (512 Fourier + Re_n + t_n) and 593,155 parameters. Train on the 51-frame LBM sequences at Re = 100 and Re = 400 with the unsteady Navier-Stokes residual (continuity + momentum with material time derivatives), an initial-condition loss pinning the rest state at t_n = 0, and the wall/lid boundary-condition loss. Optimize with 12,000 Adam steps (cosine annealing) + 1,000 L-BFGS; export the weights to pinn_temporal_re{re}.bin frame sequences and pinn_temporal_model.onnx.

Result: transient-mean u L2 of ~33% (final-frame ~30% at Re=100, ~35% at Re=400) with a peak velocity ratio of 1.13-1.16 (no overshoot); the vortex center migrates y/H ≈ 0.64 (Re=100) → 0.58 (Re=400) with Reynolds number. The surrogate animates the vortex roll-up in the cavity page's PINN Prediction panel.

Diagnosis: adding the time axis costs steady-state fidelity versus the steady-only parametric model (u L2 rises from ~24% to ~33% transient mean) -- an expected trade for continuous spatio-temporal generalization from one network. The v-field carries the larger relative error (~43-48%) because it is an order of magnitude smaller than u. This is the current project baseline and the highest-impact ML deliverable.

Iteration 5 -- Multi-scale Fourier, Re=1000, physics residuals, adaptive training

Hypothesis: the temporal baseline leaves accuracy on the table in three places -- (a) the single Fourier band (σ=5.0) mis-resolves the thin wall shear layer versus the smooth bulk, (b) the training set covers only two Reynolds numbers (100, 400), and (c) pressure is learned as a decoupled output with no PDE coupling. Each is addressable without changing the MLP core.

Change: five improvements applied together, each motivated by a prior diagnostic: (1) Multi-scale Fourier features -- three frozen sinusoid bands (σ = 1, 5, 20, m=128 each) lift (x, y) into a 768-dimensional frequency space, giving the optimizer separate purchase on the bulk recirculation and the thin wall layers. (2) Extended design space -- Re=1000 LBM data (already simulated) added to training, so the network spans laminar to transitional regimes. (3) Physics residuals -- a scale-normalized pressure-Poisson residual and vorticity-transport residual are added to the hybrid loss to couple p to the velocity structure. (4) Adaptive sensor resampling -- sensors are re-weighted by current residual every 2,000 epochs to track the moving vortex core. (5) Data pretraining -- a 1,000-epoch data-only phase runs before the hybrid loss; without it the network collapses to the constant zero-mean field (which simultaneously minimizes the data loss and the PDE residual, since u=const satisfies the steady Navier-Stokes equations). Training becomes 1,000-epoch pretrain + 8,000 Adam (cosine annealing) + 1,000 L-BFGS, 223 minutes on M5; 856,067 parameters.

Result: transient-mean u L2 of 34% (Re=100), 28% (Re=400), 33% (Re=1000); v L2 45% / 31% / 40%. The multi-scale Fourier cut the v error by ~27% relative to the single-band baseline at Re=400. The velocity standard deviation is now captured to 99% (Re=100: σ_pred = 0.0205 vs σ_true = 0.0207), versus 83% before. Final-frame u L2 is 25% (Re=100), 30% (Re=400), 38% (Re=1000).

Diagnosis: the upgrades moved the velocity field decisively, but three gaps remain. Pressure is still near-constant (pred std 0.0013 vs true 0.0413): the pressure-Poisson residual alone is insufficient because a near-constant p also nearly satisfies the Poisson equation in the bulk, so the loss gives no gradient toward the true variation -- p needs either a stronger coupling or a dedicated Poisson solver head. The Re=1000 transient is harder (final-frame u L2 38%) because its thinner boundary layer needs finer resolution than the 96x96 export grid. And the early transient (frames 0-10) is still the worst region (~47-50% L2) despite curriculum weighting, because the flow changes fastest there. The multi-scale Fourier + extended Re training is the new deployed baseline.

† Iteration 4 u L2 is the transient-mean over the 51-frame LBM sequence (final-frame ~30% at Re=100, ~35% at Re=400). Adding the time axis trades some steady-state fidelity for a single network that predicts the full evolution; see the cavity page for frame-by-frame metrics.

‡ Iteration 5 columns are transient-mean / final-frame u L2; a Re=1000 column joins the table as a new trained operating point (transient-mean u L2 33%, final-frame 38%). &Sect; Predicted pressure std remains ~0.0013 versus the true 0.0413; the pressure-Poisson residual alone was insufficient to recover the variation, as a near-constant p also approximately satisfies the Poisson equation in the bulk.

5. Results & Validation

The 3-panel comparison below shows the C++ LBM baseline, the trained PINN surrogate, and the absolute error delta for the cylinder steady-state case. The wake structure is recovered in the correct spatial location; the remaining error is concentrated in the high-gradient shear layers, exactly where the Fourier-feature work targets improvement.

The parametric cavity model demonstrates the design-space payoff directly: dragging the Reynolds-number slider from 100 to 400 shifts the primary vortex center from roughly y/H = 0.70 to y/H = 0.68, matching the known migration of the cavity recirculation with Reynolds number -- without re-running the solver. The time-parametric model goes further: a single trained network (pinn_temporal_model.onnx) animates the full vortex roll-up on the cavity page's PINN Prediction panel, interpolating continuously across both Re and time.

6. Deployment Path

The surrogate is exported with torch.onnx.export to pinn_temporal_model.onnx (and the 51-frame sequence to pinn_temporal_re{re}.bin via export_temporal.py), then served through ONNX Runtime Web so the trained field is evaluated live in the browser. This closes the loop from the C++ solver: a multi-hour simulation becomes an interactive, real-time surrogate a recruiter or engineer can probe by dragging a parameter slider or scrubbing through time. Phases 6.6 (ONNX + WASM) and 6.7 (ablation study) are complete; the temporal surrogate is the current deployed baseline.

7. Hardware & Environment

Training uses torch.device("mps") on Apple Silicon (tested on M5) -- Metal Performance Shaders, no NVIDIA GPU required. Dependencies are pinned in requirements.txt (torch, numpy, matplotlib, scipy, onnx, onnxruntime). The C++ solver is untouched: every file in this suite lives under pinn/.