Core Concepts¶

This page explains the mathematical and algorithmic foundations behind Prandtl. Each concept is traced back to its original academic source — so you understand not just how to use Prandtl, but why each component works.

1. What is a Surrogate Model?¶

A surrogate model (also called a metamodel or response surface) is a computationally cheap function $\hat{f}$ that approximates an expensive function $f$:

$$\hat{f}(\mathbf{x}) \approx f(\mathbf{x}), \quad \mathbf{x} \in \Omega \subset \mathbb{R}^d$$

In CFD, $f$ is typically a simulation that takes minutes to hours per run. The surrogate $\hat{f}$ is trained on a small set of $N$ simulation results ${(\mathbf{x}i, f(\mathbf{x}_i))}{i=1}^N$ and can predict new points in milliseconds.

Theoretical foundation: The formal framework was established by Sacks et al. (1989) in "Design and Analysis of Computer Experiments" (DACE)[^1], which treated deterministic computer simulations as realizations of a stochastic process — laying the groundwork for Gaussian Process-based surrogates. The broader field of response surface methodology (RSM) dates back to Box & Wilson (1951)[^2].

Key insight: Unlike physical experiments, computer simulations are deterministic — same inputs always produce the same outputs. This means surrogates for CFD should interpolate the training data exactly (or nearly so), rather than smoothing over noise.

[^1]: Sacks, J., Welch, W. J., Mitchell, T. J., & Wynn, H. P. (1989). Design and analysis of computer experiments. Statistical Science, 4(4), 409–423. [^2]: Box, G. E. P., & Wilson, K. B. (1951). On the experimental attainment of optimum conditions. Journal of the Royal Statistical Society: Series B, 13(1), 1–45.

2. The Four Backends¶

Prandtl provides four surrogate model backends. Each has a distinct mathematical foundation — choosing the right one depends on your data and deployment needs.

2.1 Gaussian Process (GP)¶

Theory: A Gaussian Process defines a distribution over functions. Any finite collection of function values ${f(\mathbf{x}_1), \ldots, f(\mathbf{x}_N)}$ is jointly Gaussian[^3]:

$$f \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$$

where $m(\mathbf{x})$ is the mean function and $k(\mathbf{x}, \mathbf{x}')$ is the covariance kernel. Given training data, the posterior predictive distribution at a new point $\mathbf{x}_*$ is:

$$f_ \mid \mathbf{X}, \mathbf{y}, \mathbf{x}_ \sim \mathcal{N}(\bar{f}_, \text{Var}[f_])$$

$$\bar{f}_ = \mathbf{k}_^T (K + \sigma_n^2 I)^{-1} \mathbf{y}$$

$$\text{Var}[f_] = k(\mathbf{x}_, \mathbf{x}) - \mathbf{k}_^T (K + \sigma_n^2 I)^{-1} \mathbf{k}*$$

Why this matters for CFD:

Analytical uncertainty: $\text{Var}[f_*]$ is computed exactly from the covariance structure — no sampling or approximations needed. This is what powers predict_with_uncertainty().
Interpolation: For deterministic CFD data (no noise, $\sigma_n \to 0$), the GP posterior mean passes exactly through training points.
Kernel choice encodes smoothness: The Matérn kernel family controls differentiability. Matérn-5/2 ($\nu = 2.5$) assumes twice-differentiable functions — appropriate for smooth aerodynamic coefficients.

When to use GP: Small datasets ($N < 1000$), when you need uncertainty estimates, or when you want a "safe" interpolating model.

When not to use GP: The $O(N^3)$ matrix inversion in the predictive equations makes GP impractical for $N > 5000$. GP models also cannot be exported to ONNX (non-parametric).

[^3]: Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.

2.2 Multi-Layer Perceptron (MLP)¶

Theory: An MLP with at least one hidden layer can approximate any continuous function on a compact domain to arbitrary accuracy — the Universal Approximation Theorem[^4][^5]:

For any continuous $f: [0,1]^d \to \mathbb{R}$ and $\epsilon > 0$, there exists a feedforward network with one hidden layer such that $\sup_{\mathbf{x}} |f(\mathbf{x}) - \hat{f}(\mathbf{x})| < \epsilon$.

Prandtl's MLP backend uses a shallow architecture with ReLU activations, trained with MSE loss:

$$\mathcal{L} = \frac{1}{N} \sum_{i=1}^N (\hat{f}(\mathbf{x}_i) - f(\mathbf{x}_i))^2$$

Why this matters for CFD:

Scalability: MLP training scales $O(N \cdot \text{params})$ per epoch, not $O(N^3)$. You can train on tens of thousands of points.
ONNX export: The entire model compresses to a .onnx file — deploy on edge devices, real-time control loops, or cloud inference.
GPU training: Set device='cuda' and the entire training loop runs on GPU.

Limitation: Standard MLPs provide no native uncertainty. You get point predictions only. For uncertainty with MLP, consider ensemble methods or MC Dropout[^6].

When to use MLP: Large datasets ($N > 1000$), ONNX deployment, GPU training.

[^4]: Cybenko, G. (1989). Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems, 2(4), 303–314. [^5]: Hornik, K. (1991). Approximation capabilities of multilayer feedforward networks. Neural Networks, 4(2), 251–257. [^6]: Gal, Y., & Ghahramani, Z. (2016). Dropout as a Bayesian approximation. ICML.

2.3 Random Forest (RF)¶

Theory: A Random Forest[^7] is an ensemble of $B$ decision trees, each trained on a bootstrap sample of the data with random feature subsets at each split. The prediction is:

$$\hat{f}(\mathbf{x}) = \frac{1}{B} \sum_{b=1}^B T_b(\mathbf{x})$$

Uncertainty from tree variance: Since each tree $T_b$ makes a slightly different prediction, the ensemble variance provides a natural uncertainty estimate:

$$\text{Var}[\hat{f}(\mathbf{x})] = \frac{1}{B-1} \sum_{b=1}^B (T_b(\mathbf{x}) - \hat{f}(\mathbf{x}))^2$$

Unlike GP's analytical variance, this is an empirical estimate — but it's well-calibrated for capturing epistemic uncertainty (where the model lacks data).

Key advantage: RF requires no PyTorch or GPU — pure scikit-learn, pure CPU. Train anywhere.

When to use RF: No-GPU environments, when you want tree-based uncertainty, or as a baseline before moving to GP/MLP.

[^7]: Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32.

2.4 Gradient Boosting (GB)¶

Theory: Gradient Boosting[^8] builds an additive model sequentially — each new tree corrects the residuals of the previous ensemble:

$$F_0(\mathbf{x}) = \text{argmin}\gamma \sum{i=1}^N L(y_i, \gamma)$$

$$F_m(\mathbf{x}) = F_{m-1}(\mathbf{x}) + \nu \cdot h_m(\mathbf{x})$$

where $h_m$ is a tree fitted to the negative gradient of the loss, and $\nu$ is the learning rate.

Why GB differs from RF: RF trees are independent (bootstrap + voting), while GB trees are dependent — each one learns from the mistakes of its predecessors. This makes GB more powerful but also more prone to overfitting if $\nu$ is too high or trees too many.

When to use GB: When prediction accuracy matters more than uncertainty, or as an alternative to MLP for large datasets without ONNX requirements.

[^8]: Friedman, J. H. (2001). Greedy function approximation: A gradient boosting machine. Annals of Statistics, 29(5), 1189–1232.

Decision Guide¶

Criterion	GP	MLP	RF	GB
Dataset size	< 1000	Any	Any	Any
Uncertainty	✓ Analytical	✗	✓ Ensemble variance	✗
ONNX export	✗	✓	✗	✗
GPU training	✗	✓	✗	✗
PyTorch required	✓ (GPyTorch)	✓	✗	✗
Interpolation	Exact	Approximate	Approximate	Approximate

3. Active Learning¶

Why Random Sampling is Inefficient¶

If you run 100 CFD simulations at random points in the parameter space, some regions will be oversampled while others remain sparse. The surrogate's uncertainty is uneven — tight in some areas, loose in others.

The Active Learning Principle¶

Active learning (also called sequential design in the DACE literature) asks: "Given my current surrogate, where should I run the next CFD simulation to improve it the most?"

The theory is rooted in Bayesian optimization and optimal experimental design[^9]. The key is an acquisition function $\alpha(\mathbf{x})$ that scores each candidate point. Prandtl implements the simplest and most widely-used strategy:

Maximum Standard Deviation (MaxStd):

$$\mathbf{x}{\text{next}} = \underset{\mathbf{x} \in \mathcal{X}{\text{pool}}}{\arg\max} \; \sigma(\mathbf{x})$$

where $\sigma(\mathbf{x})$ is the predictive standard deviation from the GP (or ensemble std for RF). By sampling where the model is least certain, we efficiently reduce the maximum possible error.

Alternative strategies (planned for v0.6+):

Expected Improvement (EI)[^10]: Balances exploration (uncertainty) with exploitation (predicted value). $EI(\mathbf{x}) = \mathbb{E}[\max(f_{\min} - f(\mathbf{x}), 0)]$.
Upper Confidence Bound (UCB)[^11]: $\alpha(\mathbf{x}) = \mu(\mathbf{x}) + \kappa \sigma(\mathbf{x})$, with tunable $\kappa$ controlling the exploration-exploitation trade-off.

[^9]: Chaloner, K., & Verdinelli, I. (1995). Bayesian experimental design: A review. Statistical Science, 10(3), 273–304. [^10]: Jones, D. R., Schonlau, M., & Welch, W. J. (1998). Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4), 455–492. [^11]: Srinivas, N., Krause, A., Kakade, S. M., & Seeger, M. (2010). Gaussian process optimization in the bandit setting. ICML.

4. Physics Constraints¶

CFD outputs obey physical laws — and a pure data-driven surrogate can violate them, especially in regions with sparse training data. Prandtl lets you inject domain knowledge directly into training.

4.1 Monotonicity¶

Many aerodynamic relationships are monotonic: lift coefficient CL increases with angle of attack $\alpha$ in the linear regime, thrust increases with RPM.

Theory: Monotonicity-constrained Gaussian Processes were introduced by Riihimäki & Vehtari (2010)[^12], who added virtual derivative observations:

$$P\left(\left.\frac{\partial f}{\partial x_j}\right|_{\mathbf{x}_i} \geq 0\right) \approx \Phi\left(\frac{\mu'(\mathbf{x}_i)}{\sigma'(\mathbf{x}_i)}\right)$$

where $\Phi$ is the Gaussian CDF, and $\mu'$, $\sigma'$ are the GP's predicted derivative mean and standard deviation. This penalizes negative derivatives, forcing the surrogate to respect known monotonic trends.

from prandtl import Monotonicity
surrogate.fit(X, Y, physics=[
    Monotonicity(param_idx=0, sign=1, weight=0.1)  # CL ↑ with α
])

4.2 Convexity / Concavity¶

The drag polar $C_D = C_{D0} + k \cdot C_L^2$ is convex in CL — a quadratic relationship. The boundary layer transition Reynolds number has a concave relationship with surface roughness.

Theory: Second-order derivative constraints enforce curvature direction. For a function $f$ to be convex in $x_j$, we require $\frac{\partial^2 f}{\partial x_j^2} \geq 0$ at sampled points[^13].

from prandtl import Convexity
surrogate.fit(X, Y, physics=[
    Convexity(param_idx=0, sign=1, weight=0.05)  # convex in param 0
])

4.3 Boundary Values¶

Some function values are known exactly from physics: at $\alpha = 0°$, a symmetric airfoil has exactly $C_L = 0$. At $Re \to \infty$, skin friction coefficient follows the Prandtl-Schlichting formula.

Theory: Boundary value constraints are hard equality constraints at specific points. They are implemented as high-weight penalty terms:

$$\mathcal{L}{\text{boundary}} = w \cdot |\hat{f}(\mathbf{x}{\text{boundary}}) - y_{\text{known}}|^2$$

with $w \gg 1$ (e.g., weight=10.0) to strongly enforce the constraint.

from prandtl import BoundaryValue
surrogate.fit(X, Y, physics=[
    BoundaryValue({"alpha": 0.0}, {"CL": 0.0}, weight=10.0)
])

4.4 Sobolev Training¶

Standard training matches function values only. Sobolev training[^14] also matches derivatives:

$$\mathcal{L}{\text{Sobolev}} = \mathcal{L}{\text{MSE}}(f, \hat{f}) + \lambda \cdot \mathcal{L}_{\text{MSE}}\left(\nabla f, \nabla \hat{f}\right)$$

Why this is powerful: Matching derivatives teaches the surrogate the shape of the function, not just its value at sampled points. For example, $dC_L/d\alpha = 2\pi$ (thin airfoil theory) — enforcing this gradient everywhere dramatically improves extrapolation from a few training points.

Czarnecki et al. (2017) proved that Sobolev-trained networks generalize better with fewer samples because they're learning from richer supervision — each training point contributes $(d + 1)$ pieces of information (one value + $d$ partial derivatives).

from prandtl import GradientConstraint
surrogate.fit(X, Y, physics=[
    GradientConstraint(param_idx=0, target_gradient=2 * np.pi, weight=1.0)
])

[^12]: Riihimäki, J., & Vehtari, A. (2010). Gaussian processes with monotonicity information. AISTATS. [^13]: Delecroix, M., Simar, L., & Van Keilegom, I. (2014). Convexity constraints in nonparametric regression. Handbook of Nonparametric Statistics. [^14]: Czarnecki, W. M., Osindero, S., Jaderberg, M., Świrszcz, G., & Pascanu, R. (2017). Sobolev training for neural networks. NeurIPS.

5. Multi-Fidelity Modeling (Co-Kriging)¶

The Problem¶

High-fidelity CFD (fine mesh, unsteady RANS/LES) is expensive. Low-fidelity CFD (coarse mesh, panel method) is cheap but less accurate. You have many cheap results and few expensive ones. Can you combine them?

The Solution: Co-Kriging¶

Kennedy & O'Hagan (2000)[^15] proposed an autoregressive model:

$$Z_e(\mathbf{x}) = \rho \cdot Z_c(\mathbf{x}) + Z_d(\mathbf{x})$$

$Z_c$: GP model of the cheap (low-fidelity) data
$Z_e$: GP model of the expensive (high-fidelity) data
$Z_d$: GP model of the difference between scaled cheap and expensive
$\rho$: scaling factor learned from data

The key insight: if the cheap and expensive models are highly correlated ($\rho \approx 1$), the expensive GP can "borrow strength" from the cheap data — requiring far fewer expensive runs to achieve the same accuracy.

Prandtl implementation: CoKriging accepts two datasets (X_cheap, Y_cheap) and (X_expensive, Y_expensive), and automatically learns $\rho$ and the difference GP:

from prandtl import CoKriging
ck = CoKriging(params=["alpha"], outputs=["CL"])
ck.fit(X_cheap, Y_cheap, X_expensive, Y_expensive)
Y_pred = ck.predict(X_test)

[^15]: Kennedy, M. C., & O'Hagan, A. (2000). Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87(1), 1–13.

6. Sampling Methods¶

How you choose training points directly affects surrogate accuracy. Prandtl provides three sampling strategies, each with distinct mathematical properties.

6.1 Latin Hypercube Sampling (LHS)¶

LHS[^16] divides each parameter dimension into $N$ equal-probability intervals and samples exactly one point from each interval. The result: each row and column of the 2D projection contains exactly one point.

Why it works: LHS is a stratified sampling method. It guarantees better space-filling than pure random sampling because no region is accidentally oversampled or left empty. For $d$ parameters with $N$ samples, LHS achieves variance reduction proportional to $1/N^2$ compared to $1/N$ for pure random sampling.

When to use: Default choice. Good balance of simplicity, space-filling, and statistical guarantees.

6.2 Sobol Sequences¶

Sobol sequences[^17] are low-discrepancy (quasi-Monte Carlo) sequences. Unlike random or LHS points, Sobol points are deterministic and designed to minimize the discrepancy:

$$D_N^* = \sup_{\mathbf{x} \in [0,1]^d} \left| \frac{#{\mathbf{x}i \in [0,\mathbf{x}]}}{N} - \prod{j=1}^d x_j \right|$$

This measures how uniformly the points fill the unit hypercube. Sobol sequences achieve $D_N^* = O((\log N)^d / N)$ — asymptotically better than the $O(1/\sqrt{N})$ of random sampling.

Why $N$ must be a power of 2: Sobol sequences use base-2 digit expansion. The uniformity properties are only fully realized when $N = 2^k$.

When to use: When you need maximum uniformity, and your sample count can be a power of 2. Also useful for reproducible experiments (deterministic sequence, same seed → same points).

6.3 Uniform Random¶

Simple independent random sampling from uniform distributions. No space-filling guarantees.

When to use: Baseline comparison. Not recommended for production surrogate modeling.

Method	Space-filling	Reproducible	Constraint on $N$
LHS	★★★	✓ (seed)	None
Sobol	★★★★	✓ (deterministic)	$N = 2^k$
Uniform	★	✓ (seed)	None

[^16]: McKay, M. D., Beckman, R. J., & Conover, W. J. (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21(2), 239–245. [^17]: Sobol', I. M. (1967). On the distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics, 7(4), 86–112.

7. Validation Suite¶

A surrogate model that achieves $R^2 = 0.99$ on training data is not necessarily good — it might be overfitting. Prandtl's validation suite lets you quantify how much to trust your surrogate.

7.1 Cross-Validation¶

$K$-fold cross-validation splits data into $K$ folds. For each fold, train on $K-1$ folds and test on the held-out fold. Repeat $K$ times.

$$\text{CV score} = \frac{1}{K} \sum_{k=1}^K \text{Metric}(\text{model trained without fold }k, \text{fold }k)$$

The gap between training error and CV error reveals overfitting: if CV MAE $\gg$ training MAE, the surrogate has memorized training noise rather than learning the underlying function[^18].

scores = pr.cross_validate(surrogate, X, Y, cv=5)
print(scores["CL"]["r2_mean"])    # Cross-validated R²
print(scores["CL"]["mae_std"])    # Fold-to-fold variation

7.2 Learning Curves¶

Plot model performance against training set size. A learning curve that has not plateaued means more data would improve the surrogate. A curve with a widening train-validation gap means overfitting.

$$(N, \text{MAE}{\text{train}}(N), \text{MAE}{\text{val}}(N))$$

This diagnostic answers the practical question: "Should I run 50 more CFD simulations, or is my surrogate already saturated?"

curve = pr.learning_curve(surrogate, X, Y, sizes=[20, 40, 60, 80, 100])
# Check: does val_mae keep dropping? → add more data
# Check: is train_mae << val_mae? → overfitting, simplify model

7.3 Residual Analysis¶

Residuals $r_i = y_i - \hat{y}_i$ should be:

Normally distributed — Shapiro-Wilk test, $p > 0.05$ → residuals are Gaussian (good)
Zero mean — no systematic over/under-prediction
Homoscedastic — constant variance across the prediction range

Non-normal residuals signal that the surrogate has systematic bias — it's missing some physical effect (e.g., flow separation at high $\alpha$). Skewed residuals mean the model consistently over- or under-predicts in certain regions[^19].

res = pr.residual_analysis(Y_test, Y_pred)
print(res["CL"]["shapiro_p"])      # > 0.05 → normal ✓
print(res["CL"]["skewness"])       # close to 0 → symmetric
print(res["CL"]["max_residual_idx"])  # which point is worst?

[^18]: Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed. Springer. Chapter 7. [^19]: Draper, N. R., & Smith, H. (1998). Applied Regression Analysis, 3rd ed. Wiley. Chapter 2.

8. Why CFD Surrogate Modeling Works¶

A natural question: why should a generic ML model approximate a complex CFD simulation?

The answer lies in the smoothness of aerodynamic responses. For attached flows, lift and drag coefficients are smooth, continuous functions of flow parameters ($\alpha$, Mach, Re) and geometry parameters (camber, thickness). At each point in parameter space, the Navier-Stokes equations have a unique, deterministic solution — the mapping from parameters to forces is a well-defined function, not a stochastic process.

This smoothness property is what makes surrogates viable:

GP with Matérn-5/2 kernel encodes exactly the right smoothness prior — twice differentiable, matching the physical reality of subsonic aerodynamics.
MLP with ReLU constructs a piecewise-linear approximation that can capture any smooth function (Universal Approximation Theorem).
Physics constraints correct the surrogate where data is sparse by encoding what we already know from thin airfoil theory, potential flow, and empirical correlations.

The surrogate is not "learning CFD" — it's learning the input-output mapping, which is a much simpler function than the full flow field.

Summary¶

Concept	Core Idea	Key Reference
Surrogate model	Cheap $\hat{f}$ approximates expensive $f$	Sacks et al. (1989)
Gaussian Process	Distribution over functions, analytical uncertainty	Rasmussen & Williams (2006)
MLP	Universal function approximator, ONNX exportable	Cybenko (1989)
Random Forest	Ensemble variance as uncertainty	Breiman (2001)
Gradient Boosting	Additive model of error-correcting trees	Friedman (2001)
Active Learning	Sample where uncertainty is highest	Jones et al. (1998)
Monotonicity	Virtual derivative observations	Riihimäki & Vehtari (2010)
Sobolev Training	Match values AND derivatives	Czarnecki et al. (2017)
Co-Kriging	$Z_e = \rho Z_c + Z_d$, borrow cheap-data strength	Kennedy & O'Hagan (2000)
LHS	Stratified sampling, each row/column has one point	McKay et al. (1979)
Sobol	Low-discrepancy, deterministic, $O((\log N)^d/N)$	Sobol' (1967)