4. Physics-informed graph neural networks

A conventional graph neural network learns a mapping from graph-structured inputs to outputs purely from labeled data. A physics-informed graph neural network (PI-GNN) does something stronger: it constrains that mapping so that the predicted field satisfies the governing partial differential equations, the material's constitutive response, and the prescribed boundary conditions — even in regions where no labels exist. The central idea is to treat the PDE residual, the constitutive law, and the boundary conditions as differentiable penalty terms that are evaluated on the graph itself and added to the training loss. This chapter explains how each physical constraint is embedded into the GNN architecture, how the discrete operators are constructed on an unstructured mesh graph, and how the composite loss is assembled and balanced in practice.

From mesh to graph: the discrete substrate

In solid mechanics and transport problems, the computational domain is typically discretized by a mesh. Each mesh node becomes a graph node carrying features such as coordinates x_i, a material identifier, and known boundary values. Edges connect neighboring nodes; for a finite-volume or finite-element mesh, edges are the mesh connectivity, and each edge carries a geometric vector r_ij = x_j − x_i and a length |r_ij|. The GNN's message-passing layers then act as learned stencils: node i aggregates messages from its neighbors and updates its latent state h_i. Because the graph mirrors the mesh, spatial derivatives can be approximated directly from neighbor values, which is what makes residual evaluation possible inside the network.

A practical convention must be fixed before anything else: define the outward unit normal n on each boundary face, and adopt a consistent sign convention for fluxes (for example, traction positive in tension, heat flux positive outward). Every residual and boundary term below assumes that convention, and the same convention must be used in the SVG and in the loss code.

Embedding the PDE residual

Consider a steady diffusion-type problem, ∇·(k∇u) + f = 0, where u is the unknown field, k the conductivity, and f a source term. The GNN predicts û_i at every node. To form the residual at node i, we need a discrete Laplacian. A robust choice on an unstructured graph is the weighted sum over neighbors:

L_i = Σ_j w_ij (û_j − û_i),   w_ij = k_ij · A_ij / |r_ij|
r_i = L_i + f_i

Here A_ij is the dual-face area associated with edge ij and k_ij is the edge-averaged conductivity. The residual r_i is a scalar per node and is differentiable with respect to the network weights, so it can be squared and averaged over all interior nodes to form the PDE loss L_pde = mean(r_i²). For elasticity, the same construction is applied to the momentum balance ∇·σ + b = 0, where the stress σ is recovered from the predicted displacement gradient via the constitutive law — which brings us to the second ingredient.

Embedding constitutive laws

A constitutive law maps kinematic quantities to kinetic ones: strain to stress, temperature gradient to heat flux, pressure gradient to velocity. Rather than letting the network output stress directly and hoping it is consistent with displacement, we compute the strain from the predicted displacement field using the same graph gradient operator, then apply the constitutive relation analytically. For linear elasticity with Lamé parameters λ and μ:

ε_i = 0.5 (∇û_i + ∇û_iᵀ)
σ_i = λ tr(ε_i) I + 2μ ε_i

This "hard" embedding guarantees that the predicted stress is always physically admissible for the assumed material model, and it removes an entire class of inconsistent solutions. For nonlinear materials — plasticity, hyperelasticity, or temperature-dependent conductivity — the constitutive update is written as a differentiable function and applied node-wise; the GNN learns the field, while the constitutive law remains an exact, interpretable component of the forward pass. Material parameters can be supplied as node features, allowing a single trained model to generalize across a family of materials.

Embedding boundary conditions

Boundary conditions are handled in two complementary ways. Dirichlet conditions (prescribed displacement or temperature) are enforced hard by construction: the network predicts a correction field that vanishes on Dirichlet nodes, and the known value is added back, so the output satisfies the condition exactly. Neumann conditions (prescribed traction or flux) are enforced soft through a penalty on the boundary residual:

r_bc = σ_i · n_i − t̄_i        (traction boundary)
L_bc = mean(r_bc²) over boundary nodes

Mixing hard Dirichlet and soft Neumann enforcement is standard and numerically well behaved. The boundary loss is evaluated only on tagged boundary nodes, which is why the graph must carry a boundary-type feature per node.

Assembling and balancing the composite loss

The total objective combines data fidelity (where labels exist), the PDE residual, the boundary residual, and optionally a constitutive consistency term:

L = λ_data · L_data + λ_pde · L_pde + λ_bc · L_bc

The weights λ are not cosmetic. If λ_pde dominates early in training, the network collapses toward the trivial solution û = 0, which satisfies the homogeneous PDE but violates the boundary data. A reliable strategy is to start with λ_data and λ_bc large, let the network fit the boundary and any available labels, then ramp λ_pde up over the first few thousand steps. Monitoring the individual loss components separately — not just the total — is essential: a decreasing total with a stagnant L_bc usually signals that the boundary penalty is being outvoted.

What the architecture looks like end to end

The figure below traces a single forward pass: mesh graph in, message passing, field prediction, analytic constitutive update, discrete residual and boundary evaluation, and the composite loss that back-propagates into the GNN weights.

A reproducible training example

To make the procedure concrete, consider a 2D steady heat conduction problem on a 40 × 40 structured mesh (1,600 nodes) with conductivity k = 1.0, a uniform source f = 1.0, Dirichlet u = 0 on the left edge, and Neumann flux q̄ = 0 on the other three edges. These are illustrative parameters chosen for a teaching run, not measured results. The graph is built from the mesh connectivity; each node carries (x, y, k, bc_type). A three-layer message-passing network with 64 hidden channels predicts û. The training loop evaluates the three loss terms on the graph, normalizes each by its running standard deviation, and applies the ramp schedule:

for step in range(8000):
    u_hat = model(graph)                 # hard Dirichlet applied
    L_data = mse(u_hat[labels], u_true[labels])
    L_pde  = mean(residual(graph, u_hat)**2)
    L_bc   = mean(bc_residual(graph, u_hat)**2)
    lam_pde = ramp(step, start=0.0, end=1.0, warmup=3000)
    loss = L_data + 1.0*L_bc + lam_pde*L_pde
    loss.backward(); opt.step()

What to observe: L_bc should drop quickly in the first few hundred steps as the network fits the boundary data; L_pde should fall steadily once the ramp begins; and the maximum nodal residual on interior nodes should approach the discretization error floor of the mesh (for this stencil, roughly O(h²)). A useful verification is to compare the predicted centerline temperature against a finite-volume reference solution on the same mesh — agreement within a few percent indicates the physics constraints are being honored, while a large mismatch with a small L_pde usually points to a sign error in the flux convention or a mislabeled boundary node.

The essential design principle is separation of concerns: the GNN learns the field, the constitutive law is applied analytically, and the PDE and boundary residuals are evaluated with discrete operators that mirror the mesh. Each physical constraint enters the architecture at a well-defined point, which makes the model interpretable, easier to debug, and far more data-efficient than a purely black-box surrogate.