50.4 Building a Physics Discovery Pipeline

1. The Problem: Heat Sink Design

Every smartphone, laptop, and data center server generates waste heat that, left unchecked, would destroy the very processor producing it. The small metal fins carrying that heat away conceal a surprisingly rich optimization problem. The design variables for a heat sink are the fin height \(H\), fin spacing \(S\), fin thickness \(t\), base thickness \(b\), and air velocity \(V\). The goal is to maximize total heat dissipation \(Q\) while keeping the heat sink volume below a constraint. The governing physics involves conjugate heat transfer: conduction through the metal fins coupled with convective heat transfer to the air.

Conjugate heat transfer simultaneously solves heat conduction within a solid body and convective heat transfer in the adjacent fluid. The two domains couple through shared temperature and heat flux at the solid/fluid interface. It matters because neither temperature field can be computed independently: the fin surface temperature depends on how effectively the air carries heat away, while the air temperature profile depends on how much heat the fin surface supplies. The solver iterates the energy equation in both domains until the interface temperature and heat flux converge to consistent values. Use conjugate analysis whenever the solid's thermal resistance is comparable to the convective resistance; if the solid is so conductive that its temperature is nearly uniform (Biot number, the ratio of internal conductive resistance to external convective resistance, well below 0.1), a simpler lumped or decoupled convection model suffices.

When engineers skip structured pipeline design and rely on ad hoc simulation sweeps, the result is predictable: weeks of compute exploring a fraction of the design space, with no physical insight into why one configuration outperforms another.

The full conjugate heat transfer simulation (solving coupled energy equations in the solid and fluid domains) takes 10 to 30 minutes per configuration on a workstation. Exploring the 5-dimensional design space with 10,000 evaluations would require months of compute. Our pipeline replaces this brute-force approach with four stages that complete in hours. In short: a pipeline that structures discovery into four stages replaces months of brute force with hours of principled search. The architecture of that pipeline is shown in Figure 50.4.1.

2. Stage 1: Law Discovery with PySR

We begin by running a modest number (200) of conjugate heat transfer simulations spanning the design space, and use PySR to discover the functional relationship between design variables and heat dissipation.

import numpy as np
from pysr import PySRRegressor

# Simulated dataset: 200 CFD runs varying 5 design parameters
rng = np.random.default_rng(42)
n_samples = 200

H = rng.uniform(0.01, 0.05, n_samples)      # fin height (m)
S = rng.uniform(0.002, 0.01, n_samples)      # fin spacing (m)
t = rng.uniform(0.0005, 0.003, n_samples)    # fin thickness (m)
b = rng.uniform(0.002, 0.008, n_samples)     # base thickness (m)
V = rng.uniform(0.5, 5.0, n_samples)         # air velocity (m/s)

# Physics-based synthetic data (simplified conjugate model)
# Nusselt correlation: Nu = C * Re^m * Pr^(1/3)
# with Re = V * S / nu_air
nu_air = 1.5e-5         # kinematic viscosity (m^2/s)
k_air = 0.026           # thermal conductivity (W/m/K)
k_metal = 200.0         # aluminum conductivity (W/m/K)
Pr = 0.71               # Prandtl number for air
dT = 50.0               # temperature difference (K)

Re = V * S / nu_air
Nu = 0.664 * Re**0.5 * Pr**(1.0/3.0)     # laminar flat plate
h_conv = Nu * k_air / S                    # convection coefficient

# Fin efficiency (hyperbolic tangent model)
m_fin = np.sqrt(2 * h_conv / (k_metal * t))
eta_fin = np.tanh(m_fin * H) / (m_fin * H)

# Number of fins per unit width
n_fins = 1.0 / (S + t)

# Total heat dissipation per unit depth
Q = n_fins * (2 * eta_fin * H + t) * h_conv * dT
Q_noisy = Q * (1.0 + 0.02 * rng.standard_normal(n_samples))

X = np.column_stack([H, S, t, b, V])

# Discover the law
sr_model = PySRRegressor(
    niterations=80,
    binary_operators=["+", "-", "*", "/"],
    unary_operators=["sqrt", "square", "tanh"],
    maxsize=30,
    parsimony=0.003,
    populations=40,
    population_size=60,
    X_units=["m", "m", "m", "m", "m/s"],
    y_units="W/m",
)
sr_model.fit(X, Q_noisy)

# Inspect the Pareto front
print(sr_model.equations_[['complexity', 'loss', 'equation', 'score']])

# Best equation
best_eq = sr_model.sympy()
print(f"\nDiscovered law: Q = {best_eq}")

PySR typically discovers an expression structurally similar to the analytical model: \(Q \propto \sqrt{V} \cdot H / \sqrt{S}\), capturing the key scaling relationships. The discovered expression may differ in constant factors from the full analytical model, but the functional dependencies on the design variables are correct. This is the value of symbolic regression for physics: it extracts the scaling law from data, providing physical insight that no neural network can match. PySR reports the result as a Pareto front, the set of equations that are not dominated on both complexity and accuracy, so each point offers the best accuracy achievable at its complexity level.

3. Stage 2: Parameter Calibration with PINN

The discovered symbolic law contains unknown coefficients (the constant \(C\) in the Nusselt correlation, the effective thermal conductivity). We calibrate these using experimental temperature measurements on a physical prototype. This is an inverse problem: given 20 thermocouple readings on the heat sink surface, recover the effective convection coefficient \(h_{\text{eff}}\).

import deepxde as dde
import numpy as np

# Domain: heat sink fin (2D cross-section)
# x in [0, t/2] (half-fin by symmetry), y in [0, H]
fin_thickness_half = 0.001   # 1 mm half-thickness
fin_height = 0.03            # 30 mm

geom = dde.geometry.Rectangle(
    [0, 0], [fin_thickness_half, fin_height]
)

# Unknown parameter: effective convection coefficient
h_eff = dde.Variable(50.0)  # initial guess (W/m^2/K)

# Known parameters
k = 200.0        # aluminum conductivity (W/m/K)
T_base = 80.0    # base temperature (C)
T_air = 30.0     # ambient temperature (C)

def fin_equation(x, T):
    """
    Steady-state heat conduction in a fin:
    k * (d2T/dx2 + d2T/dy2) = 0  (interior)
    """
    dT_xx = dde.grad.hessian(T, x, i=0, j=0)
    dT_yy = dde.grad.hessian(T, x, i=1, j=1)
    return k * (dT_xx + dT_yy)

# Boundary conditions
def boundary_base(x, on_boundary):
    """Base of fin: T = T_base"""
    return on_boundary and np.isclose(x[1], 0)

def boundary_tip(x, on_boundary):
    """Tip of fin: insulated (dT/dy = 0)"""
    return on_boundary and np.isclose(x[1], fin_height)

def boundary_surface(x, on_boundary):
    """Fin surface: convective BC, -k dT/dx = h(T - T_air)"""
    return on_boundary and np.isclose(x[0], fin_thickness_half)

def boundary_symmetry(x, on_boundary):
    """Symmetry axis: dT/dx = 0"""
    return on_boundary and np.isclose(x[0], 0)

bc_base = dde.icbc.DirichletBC(
    geom, lambda x: T_base, boundary_base
)
bc_symmetry = dde.icbc.NeumannBC(
    geom, lambda x: 0, boundary_symmetry
)
bc_tip = dde.icbc.NeumannBC(
    geom, lambda x: 0, boundary_tip
)

# Convective BC: -k * dT/dn = h_eff * (T - T_air)
# Implemented as a Robin BC (a boundary condition that combines
# the solution value and its normal derivative) via OperatorBC
def convective_bc(x, T, X):
    dT_dx = dde.grad.jacobian(T, X, i=0, j=0)
    return -k * dT_dx - h_eff * (T - T_air)

bc_conv = dde.icbc.OperatorBC(geom, convective_bc, boundary_surface)

# Experimental thermocouple data (20 measurements)
# Simulated from ground truth h_eff = 75 W/m^2/K
rng = np.random.default_rng(123)
observe_x = rng.uniform(0, fin_thickness_half, (20, 1))
observe_y = rng.uniform(0, fin_height, (20, 1))
observe_pts = np.hstack([observe_x, observe_y])

# Ground truth temperature field (analytical fin solution)
h_true = 75.0
m = np.sqrt(2 * h_true / (k * 2 * fin_thickness_half))
T_true = T_air + (T_base - T_air) * (
    np.cosh(m * (fin_height - observe_y))
    / np.cosh(m * fin_height)
)
T_observed = T_true + 0.5 * rng.standard_normal((20, 1))

observe_bc = dde.icbc.PointSetBC(observe_pts, T_observed, component=0)

# Assemble and solve
data = dde.data.PDE(
    geom, fin_equation,
    [bc_base, bc_symmetry, bc_tip, bc_conv, observe_bc],
    num_domain=2000,
    num_boundary=400,
    num_test=500,
)

net = dde.nn.FNN([2, 64, 64, 64, 1], "tanh", "Glorot normal")

model = dde.Model(data, net)
model.compile("adam", lr=1e-3,
              external_trainable_variables=[h_eff])
model.train(iterations=15000, display_every=3000)

model.compile("L-BFGS", external_trainable_variables=[h_eff])
model.train()

h_recovered = float(h_eff.numpy())
print(f"True h_eff:      {h_true:.1f} W/m^2/K")
print(f"Recovered h_eff: {h_recovered:.1f} W/m^2/K")
print(f"Relative error:  {abs(h_recovered - h_true)/h_true:.1%}")

In this synthetic benchmark, the PINN recovers \(h_{\text{eff}}\) with typical relative error below 3%, even with measurement noise of \(\pm 0.5^\circ\)C; accuracy on real experimental data depends on sensor placement, noise characteristics, and how well the PDE model represents the true physics. In practice, recovery error increases when the assumed PDE form omits relevant phenomena such as radiation or contact resistance. The recovered coefficient is then substituted into the symbolic law from Stage 1, producing a fully calibrated analytical model for heat sink performance. The calibrated model also informs Stage 3: the recovered \(h_{\text{eff}}\) sets the convective boundary condition in the CFD simulations that generate the FNO's training data, so the surrogate learns from physically realistic temperature fields rather than fields produced with a guessed coefficient.

4. Stage 3: Neural Operator Surrogate

The calibrated symbolic model provides the scaling law, but for detailed design optimization we need the full temperature field, not just a scalar heat dissipation value. We train an FNO surrogate on 500 conjugate heat transfer simulations, mapping design parameters to the 2D temperature distribution.

import torch
import torch.nn as nn

# Assume we have generated 500 simulation results:
# - design_params: (500, 5) design variables [H, S, t, b, V]
# - temp_fields: (500, 64, 64) temperature fields on uniform grid

def build_fno_input(design_params, resolution=64):
    """
    Construct FNO input: spatial coordinates + design parameters
    broadcast to every grid point.
    """
    n_samples = design_params.shape[0]

    # Coordinate grid
    x = torch.linspace(0, 1, resolution)
    y = torch.linspace(0, 1, resolution)
    gx, gy = torch.meshgrid(x, y, indexing='ij')

    # Broadcast design params to every grid point
    # Shape: (n_samples, 64, 64, 5 + 2)
    gx = gx.unsqueeze(0).expand(n_samples, -1, -1)
    gy = gy.unsqueeze(0).expand(n_samples, -1, -1)

    params_grid = design_params.unsqueeze(1).unsqueeze(1).expand(
        -1, resolution, resolution, -1
    )

    return torch.cat([
        gx.unsqueeze(-1),
        gy.unsqueeze(-1),
        params_grid,
    ], dim=-1)  # (n_samples, 64, 64, 7)


class HeatSinkFNO(nn.Module):
    """FNO surrogate for heat sink temperature fields."""

    def __init__(self, modes=12, width=32, n_layers=4):
        super().__init__()
        self.lift = nn.Linear(7, width)   # 2 coords + 5 params

        self.spectral_layers = nn.ModuleList()
        self.linear_layers = nn.ModuleList()
        for _ in range(n_layers):
            # Use the SpectralConv2d from Section 50.3
            self.spectral_layers.append(
                SpectralConv2d(width, width, modes, modes)
            )
            self.linear_layers.append(nn.Conv2d(width, width, 1))

        self.project = nn.Sequential(
            nn.Linear(width, 64),
            nn.GELU(),
            nn.Linear(64, 1),
        )

    def forward(self, x):
        h = self.lift(x).permute(0, 3, 1, 2)
        for spec, lin in zip(self.spectral_layers, self.linear_layers):
            h = torch.nn.functional.gelu(spec(h) + lin(h))
        h = h.permute(0, 2, 3, 1)
        return self.project(h).squeeze(-1)


# Training
surrogate = HeatSinkFNO(modes=12, width=32, n_layers=4)
optimizer = torch.optim.Adam(surrogate.parameters(), lr=1e-3)
scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, 200)

for epoch in range(200):
    surrogate.train()
    optimizer.zero_grad()

    # Forward pass (batched)
    T_pred = surrogate(fno_input_train)
    loss = torch.mean(
        torch.norm(T_pred - temp_fields_train, dim=(-2, -1))
        / torch.norm(temp_fields_train, dim=(-2, -1))
    )
    loss.backward()
    optimizer.step()
    scheduler.step()

    if (epoch + 1) % 50 == 0:
        surrogate.eval()
        with torch.no_grad():
            T_test_pred = surrogate(fno_input_test)
            test_err = torch.mean(
                torch.norm(T_test_pred - temp_fields_test, dim=(-2, -1))
                / torch.norm(temp_fields_test, dim=(-2, -1))
            )
        print(f"Epoch {epoch+1} | Train: {loss:.4f} | Test: {test_err:.4f}")

After 200 epochs, the FNO achieves a relative L2 error of approximately 2% on held-out test configurations. Each evaluation takes 5 milliseconds, compared to 15 minutes for the full computational fluid dynamics (CFD) solver: a speedup of 180,000x. The SpectralConv2d layer referenced in the model is the Fourier-space convolution defined in Section 50.3; it multiplies activations by learnable weights in the frequency domain, giving the FNO its resolution-invariant property.

5. Stage 4: Bayesian Optimization with BoTorch

With a fast surrogate in hand, we optimize the heat sink design using BoTorch, a Bayesian optimization library built on GPyTorch. Bayesian optimization is the right tool here because our surrogate, while fast, still carries approximation uncertainty. BoTorch's acquisition functions, which score how promising each candidate point is by combining the model's predicted value with its uncertainty, balance exploration (trying designs where the surrogate is uncertain) with exploitation (refining designs where the surrogate predicts high performance).

import torch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.acquisition import ExpectedImprovement
from botorch.optim import optimize_acqf
from gpytorch.mlls import ExactMarginalLogLikelihood

def objective(design_params):
    """
    Evaluate heat dissipation for a design.
    Uses the FNO surrogate for fast evaluation,
    then extracts total heat flux from the temperature field.
    """
    with torch.no_grad():
        fno_input = build_fno_input(
            design_params.unsqueeze(0), resolution=64
        )
        T_field = surrogate(fno_input).squeeze(0)

    # Total heat dissipation = integral of heat flux at base
    # Q ~ k * dT/dy|_{y=0} integrated over x
    dT_dy_base = (T_field[:, 1] - T_field[:, 0]) / (1.0 / 64)
    Q_total = torch.mean(torch.abs(dT_dy_base))

    return Q_total.unsqueeze(-1)


# Design space bounds (normalized to [0, 1])
bounds = torch.tensor([
    [0.0, 0.0, 0.0, 0.0, 0.0],   # lower bounds
    [1.0, 1.0, 1.0, 1.0, 1.0],   # upper bounds
])

# Initial sample using Sobol sequences (quasi-random points that
# fill the design space more uniformly than pseudorandom sampling)
from botorch.utils.sampling import draw_sobol_samples

n_initial = 20
X_init = draw_sobol_samples(bounds=bounds, n=n_initial, q=1).squeeze(1)
Y_init = torch.cat([objective(x) for x in X_init]).unsqueeze(-1)

# Bayesian optimization loop
n_iterations = 50
X_bo = X_init.clone()
Y_bo = Y_init.clone()

for i in range(n_iterations):
    # Fit Gaussian Process surrogate (on top of the FNO surrogate)
    gp = SingleTaskGP(X_bo, Y_bo)
    mll = ExactMarginalLogLikelihood(gp.likelihood, gp)
    fit_gpytorch_mll(mll)

    # Expected Improvement acquisition function
    best_f = Y_bo.max()
    ei = ExpectedImprovement(model=gp, best_f=best_f)

    # Optimize the acquisition function
    candidate, acq_value = optimize_acqf(
        acq_function=ei,
        bounds=bounds,
        q=1,
        num_restarts=10,
        raw_samples=256,
    )

    # Evaluate the candidate
    Y_new = objective(candidate.squeeze(0))
    X_bo = torch.cat([X_bo, candidate])
    Y_bo = torch.cat([Y_bo, Y_new.unsqueeze(-1)])

    if (i + 1) % 10 == 0:
        print(
            f"Iteration {i+1:3d} | "
            f"Best Q: {Y_bo.max():.4f} | "
            f"New Q: {Y_new.item():.4f}"
        )

# Extract the optimal design
best_idx = Y_bo.argmax()
best_design = X_bo[best_idx]
print(f"\nOptimal design (normalized): {best_design.numpy()}")
print(f"Maximum heat dissipation: {Y_bo.max().item():.4f}")

The two-level surrogate strategy (FNO as the physics model, GP as the acquisition model) is a common pattern in engineering optimization. The FNO provides a fast physics-based evaluation, while the GP captures the residual uncertainty and guides the search. After 50 iterations (each taking milliseconds), BoTorch in this benchmark finds designs within 2% of the true optimum, validated by running the full CFD solver on the top 5 candidates. Convergence speed varies with problem dimensionality and the smoothness of the objective landscape.

Practical Example: Validating the Pipeline End-to-End

We validate the complete pipeline by comparing its output against a brute-force grid search over 10,000 CFD simulations. The pipeline uses 200 simulations for PySR (Stage 1), 20 thermocouple measurements for the PINN (Stage 2), 500 simulations for the FNO (Stage 3), and 50 BoTorch iterations (Stage 4): a total of 700 simulations plus 20 experimental measurements. The brute-force approach requires 10,000 simulations. The pipeline finds a design with heat dissipation within 3% of the brute-force optimum, using 14x fewer simulations. The symbolic law from Stage 1 also provides physical insight (the \(\sqrt{V}\) velocity dependence, the \(H/\sqrt{S}\) geometry scaling) that the brute-force approach does not.

6. Extrapolation Comparison: Symbolic vs. Neural

The optimization loop finds strong designs within the training distribution, but real engineering problems rarely stay inside that distribution; new operating conditions, materials, or geometries push models into territory they have never seen.

The pipeline's symbolic and neural components have complementary strengths. We quantify this by evaluating both on an extrapolation task: predicting heat dissipation for air velocities 2x beyond the training range.

import numpy as np
from sklearn.metrics import mean_squared_error

# Training range: V in [0.5, 5.0] m/s
# Extrapolation range: V in [5.0, 10.0] m/s

# Generate extrapolation test data
V_extrap = np.linspace(5.0, 10.0, 100)
# Fix other variables at their midpoints
H_mid, S_mid, t_mid, b_mid = 0.03, 0.006, 0.0015, 0.005

Q_true_extrap = heat_dissipation_analytical(
    H_mid, S_mid, t_mid, b_mid, V_extrap
)

# Symbolic model prediction (from PySR)
Q_symbolic_extrap = sr_model.predict(
    np.column_stack([
        np.full(100, H_mid), np.full(100, S_mid),
        np.full(100, t_mid), np.full(100, b_mid),
        V_extrap,
    ])
)

# Neural network baseline (MLP trained on same data)
Q_neural_extrap = nn_model.predict(
    scaler.transform(np.column_stack([
        np.full(100, H_mid), np.full(100, S_mid),
        np.full(100, t_mid), np.full(100, b_mid),
        V_extrap,
    ]))
)

# FNO surrogate prediction
with torch.no_grad():
    fno_input_extrap = build_fno_input(
        torch.tensor(np.column_stack([
            np.full(100, H_mid), np.full(100, S_mid),
            np.full(100, t_mid), np.full(100, b_mid),
            V_extrap,
        ]), dtype=torch.float32)
    )
    Q_fno_extrap = extract_heat_flux(surrogate(fno_input_extrap))

# Compare
rmse_symbolic = np.sqrt(mean_squared_error(Q_true_extrap, Q_symbolic_extrap))
rmse_neural = np.sqrt(mean_squared_error(Q_true_extrap, Q_neural_extrap))
rmse_fno = np.sqrt(mean_squared_error(Q_true_extrap, Q_fno_extrap.numpy()))

print("Extrapolation RMSE (V = 5-10 m/s, trained on 0.5-5 m/s):")
print(f"  Symbolic (PySR):    {rmse_symbolic:.2f} W/m")
print(f"  Neural (MLP):       {rmse_neural:.2f} W/m")
print(f"  Neural Operator:    {rmse_fno:.2f} W/m")
print(f"  Symbolic advantage: {rmse_neural / rmse_symbolic:.1f}x over MLP")

The results consistently show three regimes:

The symbolic model trades interpolation accuracy for correct functional form (\(Q \propto \sqrt{V}\)), yielding the best extrapolation. The MLP achieves the tightest training fit but fails catastrophically outside the training range. The FNO falls between the two: its Fourier inductive bias provides partial extrapolation capability but cannot match a correct symbolic form.

7. Discovery Workbench Integration

Whether the pipeline follows today's modular pattern or tomorrow's end-to-end differentiable approach, it needs a software harness that tracks each stage's status, stores intermediate artifacts, and coordinates execution.

The four-stage pipeline integrates naturally into the Discovery Workbench architecture from Chapter 6. Each stage maps to a workbench component:

from dataclasses import dataclass, field
from typing import Optional
from enum import Enum

class StageStatus(Enum):
    PENDING = "pending"
    RUNNING = "running"
    COMPLETED = "completed"
    FAILED = "failed"

@dataclass
class PhysicsDiscoveryPipeline:
    """
    Four-stage physics discovery pipeline
    integrated with the Discovery Workbench.
    """
    # Stage 1: Law Discovery
    symbolic_model: Optional[object] = None
    discovered_law: Optional[str] = None
    pareto_front: list = field(default_factory=list)

    # Stage 2: Parameter Calibration
    calibrated_params: dict = field(default_factory=dict)
    calibration_error: float = 0.0

    # Stage 3: Surrogate Model
    surrogate_model: Optional[object] = None
    surrogate_error: float = 0.0

    # Stage 4: Optimization
    optimal_design: Optional[dict] = None
    optimal_objective: float = 0.0
    optimization_history: list = field(default_factory=list)

    # Pipeline metadata
    stage_status: dict = field(default_factory=lambda: {
        "law_discovery": StageStatus.PENDING,
        "calibration": StageStatus.PENDING,
        "surrogate": StageStatus.PENDING,
        "optimization": StageStatus.PENDING,
    })

    def run_law_discovery(self, X, y, units_X, units_y):
        """Stage 1: Discover governing law with PySR."""
        self.stage_status["law_discovery"] = StageStatus.RUNNING

        from pysr import PySRRegressor
        model = PySRRegressor(
            niterations=80,
            binary_operators=["+", "-", "*", "/"],
            unary_operators=["sqrt", "square", "tanh"],
            maxsize=30,
            parsimony=0.003,
        )
        model.fit(X, y, X_units=units_X, y_units=units_y)

        self.symbolic_model = model
        self.discovered_law = str(model.sympy())
        self.pareto_front = model.equations_.to_dict('records')
        self.stage_status["law_discovery"] = StageStatus.COMPLETED

        return self.discovered_law

    def run_calibration(self, observe_pts, observe_vals, pde_fn):
        """Stage 2: Calibrate parameters with PINN."""
        self.stage_status["calibration"] = StageStatus.RUNNING
        # ... PINN inverse problem (as in Stage 2 above) ...
        self.stage_status["calibration"] = StageStatus.COMPLETED
        return self.calibrated_params

    def run_surrogate(self, X_train, y_train, X_test, y_test):
        """Stage 3: Train FNO surrogate."""
        self.stage_status["surrogate"] = StageStatus.RUNNING
        # ... FNO training (as in Stage 3 above) ...
        self.stage_status["surrogate"] = StageStatus.COMPLETED
        return self.surrogate_error

    def run_optimization(self, bounds, n_iterations=50):
        """Stage 4: Bayesian optimization over surrogate."""
        self.stage_status["optimization"] = StageStatus.RUNNING
        # ... BoTorch optimization (as in Stage 4 above) ...
        self.stage_status["optimization"] = StageStatus.COMPLETED
        return self.optimal_design

    def report(self):
        """Generate a summary of the discovery pipeline."""
        return {
            "discovered_law": self.discovered_law,
            "calibrated_params": self.calibrated_params,
            "surrogate_error": f"{self.surrogate_error:.2%}",
            "optimal_design": self.optimal_design,
            "optimal_objective": self.optimal_objective,
            "stages_completed": sum(
                1 for s in self.stage_status.values()
                if s == StageStatus.COMPLETED
            ),
        }

# Usage
pipeline = PhysicsDiscoveryPipeline()
law = pipeline.run_law_discovery(X, Q_noisy,
    units_X=["m", "m", "m", "m", "m/s"],
    units_y="W/m",
)
print(f"Discovered: Q = {law}")
print(f"Pipeline status: {pipeline.report()}")

The pipeline object tracks the status of each stage, stores intermediate results, and generates summary reports. In a production Discovery Workbench, this pipeline would be registered as a discovery agent (see Chapter 53), capable of running autonomously and reporting results to the experiment registry (see Chapter 47).

8. Lessons and Decision Framework

With the pipeline built and its components integrated into a trackable workbench, the remaining question is practical: given a new physics or engineering problem, which of these methods should you reach for first?

Three core methods address physics and engineering discovery. Choosing among them depends on three questions:

Three Decision Questions

Do you need interpretability? If yes, start with symbolic regression (Section 50.1). The discovered equations provide physical insight, extrapolate correctly, and can be published in papers. Use PySR with dimensional constraints for any problem where you believe a compact law exists.

Do you need to recover unknown parameters? If yes, use PINNs (Section 50.2). They excel at inverse problems where the PDE form is known but coefficients are missing. Combine with Noether's theorem priors (the principle that every continuous symmetry of a physical system corresponds to a conserved quantity) for conservative systems, or simulation-based inference (SBI) for problems with intractable likelihoods.

Do you need fast evaluation across many configurations? If yes, train a neural operator surrogate (Section 50.3). FNO for regular grids, DeepONet for irregular sensor locations. The break-even point is roughly 1,000 evaluations; below that, the classical solver is cheaper.

For complex engineering design problems, use all three in sequence (this section). The symbolic model provides the scaling law and physical understanding. The PINN calibrates the model to experimental reality. The neural operator enables rapid design-space exploration. And Bayesian optimization finds the optimum efficiently. This four-stage pattern, general enough for any physics or engineering domain, is the core recipe of this chapter.