The hype and the reality: Why graph neural networks struggle with sharp physics

Credit: Wikipedia

The GNN Gold Rush

Graph Neural Networks are having a moment. From molecular dynamics to weather forecasting to traffic prediction, the pitch is always the same: Your data has structure, let a GNN exploit it. And for many problems, this pitch delivers. GNNs have produced state-of-the-art results across dozens of benchmarks precisely because they encode a powerful inductive bias which is that nearby things influence each other.

In scientific computing and physics simulation, this reasoning sounds especially compelling. Partial differential equations (PDEs), the mathematical language of physics, are inherently local. The Laplacian at a point depends only on its infinitesimal neighbourhood. Finite element methods have exploited mesh locality for decades. So when the deep learning community asked can we replace meshes with graphs and solvers with neural networks?, the answer seemed obvious.

We decided to test this intuition with a focused experiment. We took one of the simplest and most well-studied PDEs in computational physics, the 1D viscous Burgers’ equation, and systematically replaced more and more of a standard neural network solver with a Graph Neural Network. What we found was not what we (or at least I) expected.

The GNN didn’t fail because of engineering issues, hyperparameter choices, or insufficient training. It failed because the mathematical operation at the heart of all GNNs, message passing, is fundamentally at odds with the physics it was trying to capture.

This is not a new theoretical observation (the overs-moothing problem in GNNs is well-documented [1]), but seeing it play out in the context of PDE solving, where we can compare against exact solutions, makes the failure mode concrete and quantifiable.

Setting the Stage: PINNs and Burgers’ Equation

What is a PINN?

A Physics-Informed Neural Network (PINN) approximates the solution of a PDE by a neural network. Rather than training on simulation data, PINNs embed the governing equations directly into the loss function via automatic differentiation. The optimisation problem is:

Eq 1: The optimisation problem with a physics loss function

The standard architecture for PINNs is a Multi-Layer Perceptron (MLP), a fully connected network that processes each space-time point (x,t) independently. By the universal approximation theorem, a sufficiently wide single-hidden-layer network can approximate any continuous function on a compact set. In practice, moderate-depth tanh networks work remarkably well.

Why Burgers’ Equation?

The 1D viscous Burgers’ equation is a perfect test case:

Eq 2: Burgers equation in 1D

I explored the PiNN implementation for solving this in a pervious blog [2]. This equation has a beautiful and physically important feature. Writing it in conservation form:

Eq 3: Conservative form of Burgers equation

we see that the flux f(u)=u²/2, is nonlinear. Characteristics carry information at speed uu itself — regions where u>0 propagate rightward while u<0 regions propagate leftward. The wave profile steepens until adjacent characteristics cross, at which point a shock (discontinuity) forms.

In the inviscid limit (ν→0), the Rankine–Hugoniot condition gives the shock speed:

Eq 4: Shock speed in burgers equation

For asymmetric initial condition u(x,0) = −sin⁡(πx), the states on either side of the shock satisfy u_L​ = −u_R​, giving s = 0, the shock sits permanently at x = 0.

With finite viscosity ν = 0.01/π ≈ 3.18×10−3, the discontinuity is smoothed into a thin internal layer of width O(ν). But this layer is approximately 3 grid points wide on a typical computational mesh, making it a severe test for any method that involves spatial averaging.

This is why we chose Burgers’ equation: it has smooth regions where spatial correlation is useful and a sharp shock where it’s destructive. It’s the ideal litmus test for GNNs.

Experiment 1: The MLP Baseline

Before adding any graph structure, we establish a baseline. We simply present the results of my previous blog here [1] as reference where a simple neural network was used with the physics loss function.

Fig 1:The velocity profile of the fluid showing the shock cleanly at x = 0 [1]

The MLP captures the shock cleanly. The solution is visually indistinguishable from the exact Hopf–Cole analytical solution. This is our gold standard.

Experiment 2: Pure GNN. No MLP at All

With our baseline established, we went straight for the bold hypothesis: replace the MLP entirely with a GNN. If spatial structure matters for PDEs, and GNNs encode spatial structure, then a GNN-only PINN should work at least as well as an MLP, right?..

Architecture

We built a pure GNN-PINN with no MLP component whatsoever. The architecture is: a 2-layer encoder that lifts (x,t) into a hidden space, 3 GraphConv message-passing layers with residual connections, and a linear decoder:

[embed]

Note what’s missing: there is no independent pointwise pathway. Every piece of the prediction flows through the message-passing layers. The model has 17,901 parameters, more than double the MLP’s 7,701 so it’s not starved for capacity.

We construct per-time-slice KNN chain graphs (64 spatial nodes × 40 time slices = 2,560 collocation points).

The Physics Loss

The physics residual is computed through automatic differentiation:

[embed]

Training Strategy

We use the standard two-phase approach that has become convention in the PINN literature:

1. Adam (8,000 steps, lr = 10−3): stochastic gradient descent for fast initial convergence
2. L-BFGS (4,000 iterations, strong Wolfe line search): quasi-Newton method for high-precision refinement

The shift from Adam to L-BFGS is crucial. Adam navigates the loss landscape quickly but struggles with the final 10–4 to 10−7 precision. L-BFGS, which builds an approximate Hessian Bk​ ≈ ∇2L, achieves superlinear convergence in this regime.

Results

The first thing to notice is the training instability. The Adam loss oscillates wildly, going up from epoch 1000 to 2000, then back down, then up again at 8000. Compare this to the MLP, whose Adam trajectory was a smooth monotonic decrease. The GNN’s loss landscape is significantly more rugged, likely because small weight changes alter the message-passing dynamics globally rather than locally.

But L-BFGS eventually forces convergence, and the final loss of 3.5×10−7 looks excellent on paper. Let’s see the plot:

Fig 2: The velocity profile from a only GNNs

The solution is visibly diffuse. The shock, which should be a razor-thin transition at x = 0, is smeared into a broad gradient. The colour transition from blue to red is gradual rather than sharp. The smooth regions away from the shock look approximately correct, but the critical feature, the discontinuity, is destroyed. This was our first indication that something fundamental was wrong with the GNN approach for this problem. But we weren’t ready to give up yet, perhaps the GNN just needed help from an MLP to handle the sharp bits.

Experiment 3: The Hybrid Approach. MLP + GNN (GraphConv at α ≈ 30%)

The pure GNN’s failure motivated a different strategy: rather than replacing the MLP, augment it. Let the MLP handle what it’s good at (sharp features) and add a GNN branch to contribute spatial coherence. We use a hybrid architecture with a learnable mixing parameter.

The Hybrid Model

The MLP and GNN run in parallel. Their outputs are blended via a learnable parameter α∈[α_min⁡,1] :

Eq 5: Combined GNN and MLP outputs

where α = max⁡ ⁣(σ(αraw), αmin⁡) and αraw∈R is optimised jointly with the network weights.

The GNN Branch

The GNN branch uses GraphConv from PyTorch Geometric — a message-passing layer with un-normalised sum aggregation:

Eq 6: Message passing in GNNs

We stack 2 such layers with residual connections, preceded by a 2-layer encoder and followed by a linear decoder:

[embed]

Graph Construction

For 1D spatial data, we construct a k-nearest-neighbour chain graph: each spatial node connects to its k closest neighbours along x. One such graph is built per time slice, and they are batched via edge index offsetting:

[embed]

Training Protocol

We use a two-phase curriculum:

1. Phase 1: Train the pure MLP to convergence (identical to baseline)
2. Phase 2: Create the hybrid model, copy the converged MLP weights, and fine-tune with the GNN active

This avoids the GNN corrupting the MLP’s solution during early, volatile training.

Result

Fig 4: Final loss results

Fig 5: We see that the GNN + MLP sol seems much closer to the actual solution with the shock being captured well in the centre, even if its not as sharp as the ground truth result

This looks great. The shock is sharp, the physics is correct, and the loss actually improved over the pure MLP by a factor of ~30. The network chose α≈0.29, allowing the GNN to contribute about 30% of the prediction while the MLP retains primary control.

At first glance, this validates the GNN hypothesis. The spatial inductive bias genuinely helps, as a supplement, not a replacement. The hybrid model achieves a 30× lower loss than the MLP alone, and beats the pure GNN’s visual quality by a wide margin. But what happens if we ask the GNN to do more?

Experiment 4: The GNN Takes Over (GraphConv at α ≈ 70%)

We now build a GNN-dominant architecture.

Fig 6: Dominant GNN architecture

The deeper, wider GNN now carries at least 50% of the prediction, with a 4-layer message-passing backbone:

[embed]Fig 7: Schematic of how the heavy gnn is built

Results

Fig 7: Final loss and alpha values from training

The loss is lower than any previous experiment. On paper, this is our best model. Now look at the plot:

Fig 8: The velocity profile with a dominantlt Graph PiNN

The shock is destroyed.

Horizontal banding artefacts appear everywhere. The transition that should be a razor-thin front spanning O(ν)≈3 grid points is now a gradient smeared across nearly half the domain. The low loss tells us the model satisfies the PDE in a least-squares sense, but it has learned a fundamentally wrong solution — one that minimises the residual by flattening the spatial gradients, rather than resolving them.

This is the central point. Let’s understand why.

Why GNNs Fail at Shocks: A Mathematical Argument

The Smoothing Effect of Sum Aggregation

The GraphConv layer performs message passing as shown in Eq 6 above.

The second term is a weighted sum of neighbour features. In matrix form, this is multiplication by the adjacency matrix:

Eq 7: Message passing in matrix form

If we ignore the nonlinearity σ and the self-transform Θ1​ for a moment, L layers of message passing produce:

Eq 8: Linear terms of the message passing updates (this is ofc a simplification as the non-linear terms will add effects, but at least for small amplitude, in a taylor expansion, this would be the leading term)

a polynomial in the adjacency matrix applied to the initial features.

To understand what this polynomial does to a signal, we need to build up the spectral machinery step by step. The payoff is a precise understanding of why GNNs act as low-pass filters and hence supress shocks.

Step 1: The adjacency matrix as message passing

Recall that the adjacency matrix A has entries Aij=1 if nodes i and j are connected, 0 otherwise. Multiplying A by a signal vector h=(h1,…,hN)^T gives:

Eq 9: Message passing sum betweek neighbours

This is exactly the message-passing sum, node i receives the sum of its neighbours’ values. So the question “what do L layers of message passing do?” reduces to “what does multiplying by A repeatedly do to a signal on the graph?”

Step 2: The degree matrix and normalisation

The degree of node i is di=∑j Aij(its number of neighbours). The degree matrix D is diagonal with Dii=di

On our 1D chain graph with k neighbours on each side, interior nodes have degree 2k while boundary nodes have fewer. This non-uniformity is a problem: multiplying by A scales the output by the node degree, making eigenvalues depend on local connectivity rather than on the signal’s frequency content.

The fix is symmetric normalisation. Define:

Eq 10: Normalisation of adjaceny matrix

which has entries A^ij=Aij/ \sqrt(di dj​​). This ensures that the resulting matrix has eigenvalues in a controlled range, regardless of the degree distribution.

Step 3: The graph Laplacian and its eigendecomposition

The normalised graph Laplacian is:

Eq 11: Normalised graph lapacian operator

Applied to a signal, (\tilde{L}h)_i​ measures how different node i’s value is from its (normalised) neighbours. It is the discrete analogue of −∇2h on a continuous domain. If h is constant, \tilde{L}h=0. If h oscillates rapidly, \tilde{L}h is large.

Since \tilde{L} is real and symmetric, the spectral theorem guarantees a complete orthonormal eigenbasis:

Eq 12: Spectral decomposition of laplacian operator

where ϕk​ are orthonormal eigenvectors and the coefficients have the following properties

Eq 13: Spectral coefficients properties

Step 4: Eigenvectors as graph Fourier modes

The eigenvectors ϕk​ are the graph analogue of Fourier modes. On a regular 1D chain, they literally are the discrete cosine modes. The eigenvalue λ~k measures the frequency of the mode:

Fig 9: Schematic of eigenvectors and their spectral interpretation

Why λ~0=0? The constant vector 11 satisfies L1=0 because (Lh)i=∑j(hi−hj)=0 when all values are equal. A constant signal has zero variation it’s the smoothest possible mode.

The Rayleigh quotient makes the frequency interpretation precise:

Eq 14: Frequency intepretation of the modes

which sums the squared differences of the eigenvector across all edges. Small λ~k= slowly varying (smooth). Large λ~k​ = rapidly oscillating.

Step 5: The Graph Fourier Transform

Any signal can be expanded in this eigenbasis:

Eq 15: Fourier expansion in eigenbasis of the laplace operator

Since L~=I−A, the normalised adjacency has the same eigenvectors but with eigenvalues flipped:

Eq 16: Adjacency matrix in terms of the laplacian

Applying the A matrix to the signal gives:

Eq 17: Applying the Adjanceny operator to the signal

Each Fourier coefficient h^k gets multiplied by the factor μ_k=(1−λ_k):

Fig 10: Intepreting the adjacency factor with each eigenvalue

Step 7: Repeated application, multiple GNN layers

After L layers of message passing (ignoring nonlinearities), the signal has been multiplied by A^L:

Eq 18: Repeated application of A

Fig 11

The key observation: ∣μk∣=∣1−λ~k∣< for all λ~k∈(0,2). Raising to the LL-th power sends these terms to zero exponentially fast. Only the λ~0=0 mode (the constant / DC component) survives with factor 1L=1. Everything else is killed.

This is exponential suppression of all non-constant modes which is exactly what a low-pass filter is, that gets more aggressive with every layer. With L=4 layers, any mode with λ~k>0.5 has been multiplied by at most 0.0625 = 0.5⁴, a 94% reduction.

Step 8: Why this kills shocks

A shock is a near-discontinuity. Its Graph Fourier decomposition has significant energy at all frequencies (coefficients h^k∼1/k). The GNN’s low-pass filter systematically kills the high-k components that make the step sharp. After L layers, the high-frequency content needed to represent the sharp transition has been exponentially suppressed. The shock becomes a gentle slope.

Step 9: The polynomial filter view

In practice, GNNs don’t just compute A^L, they have learnable weights clcl​ at each layer, so the effective operation is a polynomial filter:

Eq 19

The GNN can learn any degree-L polynomial of A, but a degree-L polynomial can have at most L roots. With L=4 layers, the frequency response p(μ) can cross zero at most 4 times. It cannot selectively preserve an arbitrary set of high frequencies while suppressing others. The filter is inherently smooth and low-resolution in the spectral domain. To preserve a shock (which requires p≈1 for all modes), the GNN would need c0=1,c_{l>0}=0, in other words, no message passing at all.

Step 10: Why nonlinearities don’t save you

The analysis above is exact in the linear regime (no activation functions). Real GNN layers have tanh activations between them, so the computation isn’t literally A^L. However, each layer still averages over neighbours, which smooths the signal. The nonlinearity can reshape magnitudes but cannot undo the information mixing that has already occurred. The oversmoothing literature (Li et al., 2018; Oono & Suzuki, 2020) confirms that deep GNNs with nonlinearities still converge node features toward a common value and our experiments confirm it for this specific PDE.

The Numbers Don’t Lie. But They Don’t Tell the Whole Truth

Here is the paradox that makes this result so important:

Fig 12: A summary of the results from the different experiments

The two lowest-loss models are the two worst models.

The pure GNN (Experiment 2) and the GNN-dominant hybrid (Experiment 4) both achieved losses below 5×10−7, an order of magnitude better than the pure MLP. Yet both produced physically wrong solutions with smeared shocks.

How is that possible? Because the physics loss is an average over collocation points. A solution that is “roughly correct” everywhere, with gentle gradients that approximately satisfy the PDE, can achieve a low average squared residual. In fact, a smoother solution has smaller spatial derivatives, making the PDE residual easier to satisfy everywhere except at the shock. The GNN achieves a low loss precisely because it smoothed the shock. The smoothing helps the loss, even though it destroys the physics.

Consider the shock region. The true solution changes by Δu∼2over a distance Δx∼0.003, giving ∣ux∣∼600. The physics loss must balance terms of this magnitude. A smooth approximation with ∣ux∣∼5t rivially satisfies the PDE because all the troublesome nonlinear term s are small. It’s wrong, but it’s easily wrong.

This is a known issue in PINN training: the L2 loss does not penalise local errors at sharp features proportionally to their physical importance. The GNN’s smoothing makes this problem worse by systematically removing the high-frequency content that the loss should penalise but averages away.

The lesson: loss values are not sufficient metrics for PDE solvers. You must inspect the solution. If someone tells you their Graph-PINN achieved a lower loss than an MLP-PINN, ask to see the plot.

What the Network Itself Tells Us

Perhaps the most revealing finding is what happens when we let the network choose α freely (no floor constraint).

In Experiment 3, with αmin⁡=0.15, the network converged to α=0.2946. It could have pushed α higher, there was no ceiling, but it chose to keep the GNN’s contribution below 30%.

We ran this experiment multiple times with different random seeds. The network consistently settled in the range α∈[0.27,0.32] The optimiser independently discovered that ~30% GNN is the sweet spot.

When we forced α≥0.50 (Experiment 4), the loss improved but the solution degraded. The network was telling us the right answer through its learned α, and we overrode it.

This is a general principle worth remembering: when a learnable parameter consistently converges to a particular value across runs, that value is informative. The gradient of the loss with respect to ααencodes the marginal benefit/cost of shifting prediction weight between architectures. The equilibrium α∗ represents the point where:

Eq 20: The optimisation of parameters needs to be analsed from both cost and benefit

At α∗≈0.30, the spatial regularisation from the GNN (benefit) exactly balances its smoothing over the shock (cost).

The Broader Picture: When Does This Matter?

Problems Where GNNs Will Struggle

Our findings generalise beyond Burgers’ equation. Any PDE whose solutions contain sharp gradients, discontinuities, or thin layers will be adversely affected by message-passing architectures:

• Compressible fluid dynamics: shocks, contact discontinuities, expansion fans
• Reaction-diffusion systems: travelling waves with sharp fronts (flame sheets, chemical reaction fronts)
• Phase-field models: interfaces between phases (width ∝ϵ→0∝ϵ→0)
• Elastic-plastic deformation: yield surfaces and strain localisation bands
• Hyperbolic conservation laws generally: any equation admitting weak solutions with discontinuities

The common thread is the presence of features whose characteristic width is comparable to or smaller than the GNN’s receptive field. When the GNN “sees across” a discontinuity, it averages across it.

Problems Where GNNs Excel

Conversely, GNNs are well-suited to problems where the solution is smooth and spatial coherence is the dominant feature:

• Elliptic PDEs (Poisson, Laplace): solutions are C∞ in the interior
• Diffusion-dominated flows: the physics is itself a smoothing operation
• Steady-state problems: no time evolution means no shock formation
• Irregular geometries: graphs naturally represent complex domains without structured meshes
• Multi-body systems: edges encode interactions (molecular dynamics, particle systems)

The Meta-Lesson

The GNN hype cycle has produced a cognitive bias in the community: if your data has spatial structure, you should use a GNN. Our experiment provides a concrete counterexample, one that illustrates a well-understood theoretical limitation (oversmoothing) in a setting where it has direct physical consequences.

The right question is not “does my data have structure?” but “does my problem benefit from spatial averaging?” If the answer is no, if the interesting physics lives in sharp gradients rather than smooth fields, then a GNN will systematically destroy the very features you need.

A 7,701-parameter MLP with tanh⁡tanh activations captured a shock perfectly. A 17,901-parameter pure GNN smeared it. A 25,000+ parameter GNN-dominant hybrid, with carefully designed message passing, LayerNorm, and residual connections, smeared it into oblivion.

More parameters and more sophisticated architecture does not always mean better. Sometimes, the simplest tool is the right one.

Conclusion

We set out to test whether Graph Neural Networks could improve Physics-Informed Neural Networks for solving PDEs. The answer is nuanced:

As a standalone architecture, a pure GNN fails. Even with generous training (12,000 total iterations), the pure GNN produces a diffuse solution with a smeared shock. Message passing cannot represent the sharp discontinuity, period.

At low doses (~30%), GNNs provide a useful spatial inductive bias that can improve loss convergence without degrading solution quality. The GNN acts as a spatial regulariser in smooth regions, complementing the MLP’s ability to capture sharp features.

At high doses (>50%), GNNs destroy the solution. The mathematical mechanism is clear: message passing is equivalent to applying a low-pass graph filter, smoothing sharp features into gentle gradients.

The loss function is complicit. The L2 physics loss averages over collocation points, allowing smooth-but-wrong solutions to achieve low residuals. The GNN-dominant model had the lowest loss of any experiment while producing the worst solution.

The takeaway for practitioners: don’t assume GNNs are always beneficial just because your data has spatial structure. Profile the solution, not just the loss. And when someone claims their Graph-PINN outperforms an MLP-PINN, ask one simple question: “Can I see the plot?”

Darsh Kodwani is on LinkedIn

References

[1] https://arxiv.org/pdf/2006.13318

[2] https://medium.com/@darshkodwani13/a-simple-application-of-pinns-burgers-equation-3a0ba83e3904