Divergence-Free Neural Networks

• dfNNs are neural parameterizations that enforce divergence-free conditions exactly through architectural designs, ensuring conservation laws hold identically.
• They leverage methods like potential lifting and spectral projection to transform physical constraints into well-defined neural representations, enhancing stability and accuracy.
• Applications in incompressible flows, solid mechanics, and medical imaging demonstrate that encoding physical laws directly in the network improves model performance.

Divergence-Free Neural Networks (dfNNs) are neural parameterizations whose outputs satisfy a divergence-free constraint identically, by construction, rather than only approximately through a loss penalty. In the literature, the constraint appears in several equivalent or application-specific forms: the incompressibility condition for velocity fields, $\nabla \cdot u = 0$; the equilibrium condition for stress fields, $\mathrm{div}\,P = 0$; and the space-time conservation law $\partial_t \rho + \nabla \cdot (\rho v)=0$, recast as divergence-freeness of a space-time current $J=(\rho,\rho v)$ in $\mathbb{R}^{n+1}$ (Richter-Powell et al., 2022, He et al., 30 Mar 2026). The defining feature of dfNNs is therefore architectural: the admissible hypothesis space is restricted to solenoidal fields, or to divergence-free tensor fields, before optimization begins.

1. Mathematical foundations

A central starting point is the observation that the continuity equation can be rewritten as a divergence-free condition in space-time. If

$J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$

then

$$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$

Parameterizing $J$ as a divergence-free object therefore enforces conservation exactly. In the differential-forms formulation, this is achieved by identifying the metric-lowered current $J^\flat$ as a co-closed $1$-form and representing it as a coexact form, $\mathrm{div}\,P = 0$0, or equivalently as a curl-like construction in space-time (Richter-Powell et al., 2022).

The same idea admits a classical vector-calculus interpretation. In Euclidean space, for a smooth vector field $\mathrm{div}\,P = 0$1 with Jacobian $\mathrm{div}\,P = 0$2, the divergence is $\mathrm{div}\,P = 0$3. For linear fields $\mathrm{div}\,P = 0$4, skew-symmetry $\mathrm{div}\,P = 0$5 implies $\mathrm{div}\,P = 0$6, hence divergence-freeness. More generally, divergence-free vector fields are linked to skew-adjoint directional derivative operators on $\mathrm{div}\,P = 0$7 spaces, and the associated flows preserve volume by Liouville’s theorem; orthogonal recurrent models exploit precisely this structure in discrete time (Azencot et al., 2021).

For Sobolev vector fields on bounded Lipschitz domains, the modern dfNN literature makes this structure explicit through potential lifting. On domains that are starlike with respect to a ball, every divergence-free Sobolev field $\mathrm{div}\,P = 0$8 admits a stable potential $\mathrm{div}\,P = 0$9 such that

$\partial_t \rho + \nabla \cdot (\rho v)=0$0

with

$\partial_t \rho + \nabla \cdot (\rho v)=0$1

and

$\partial_t \rho + \nabla \cdot (\rho v)=0$2

In two dimensions this reduces to the stream-function representation $\partial_t \rho + \nabla \cdot (\rho v)=0$3, while in three dimensions it reduces to the vector-potential representation $\partial_t \rho + \nabla \cdot (\rho v)=0$4 (He et al., 30 Mar 2026).

2. Architectural realizations

The most direct dfNN construction parameterizes a potential and differentiates it. In two dimensions, one sets

$\partial_t \rho + \nabla \cdot (\rho v)=0$5

and in three dimensions

$\partial_t \rho + \nabla \cdot (\rho v)=0$6

Because $\partial_t \rho + \nabla \cdot (\rho v)=0$7 and $\partial_t \rho + \nabla \cdot (\rho v)=0$8, incompressibility holds identically. A higher-dimensional variant uses an antisymmetric tensor potential $\partial_t \rho + \nabla \cdot (\rho v)=0$9 and defines the velocity row-wise by

$J=(\rho,\rho v)$0

Then

$J=(\rho,\rho v)$1

in the distributional sense because $J=(\rho,\rho v)$2 and mixed derivatives commute (He et al., 30 Mar 2026, Cheng et al., 18 Mar 2026).

A second family enforces divergence-freeness by projection. On periodic domains, the Leray projector

$J=(\rho,\rho v)$3

maps any field to its divergence-free component. In Fourier variables, for $J=(\rho,\rho v)$4,

$J=(\rho,\rho v)$5

with the zero mode removed to enforce zero mean. This produces a differentiable spectral layer that can be inserted into deterministic neural operators or into generative probability flows; exact incompressibility then holds up to discretization error (Li et al., 25 Mar 2026, Genuist et al., 27 Jan 2026). A related discrete perspective appears in large-eddy simulation, where divergence-freeness is enforced not by the network alone but by a divergence-consistent filter and projection pipeline acting on the full right-hand side of the coarse dynamics (Agdestein et al., 2024).

The same architectural logic extends from vectors to tensors. In solid mechanics, a stress row $J=(\rho,\rho v)$6 can be represented as

$J=(\rho,\rho v)$7

or, in Fourier space,

$J=(\rho,\rho v)$8

which guarantees $J=(\rho,\rho v)$9 mode-by-mode (Khorrami et al., 2024, Khorrami et al., 1 May 2026). For divergence-free symmetric tensors, Riemann Tensor Neural Networks (RTNNs) represent

$\mathbb{R}^{n+1}$0

with $\mathbb{R}^{n+1}$1 satisfying Riemann symmetries; the output is symmetric and divergence-free to machine precision by construction (Jnini et al., 2 Mar 2025).

3. Approximation theory and optimization regimes

A distinctive theoretical development is the divergence-free finite neuron space built from linearized ReLU$\mathbb{R}^{n+1}$2 features. With the augmented variable $\mathbb{R}^{n+1}$3 and quasi-uniform inner parameters $\mathbb{R}^{n+1}$4, the scalar finite neuron space is

$\mathbb{R}^{n+1}$5

where $\mathbb{R}^{n+1}$6. Combining scalar Sobolev integral representation with stable antisymmetric potential lifting yields divergence-free basis functions

$\mathbb{R}^{n+1}$7

and the approximation theorem

$\mathbb{R}^{n+1}$8

In particular, the $\mathbb{R}^{n+1}$9 rate is

$J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$0

and the $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$1 rate is available only for $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$2 (He et al., 30 Mar 2026).

Optimization strategies vary with the architecture. In finite-neuron and TransNet-style constructions, inner parameters are fixed and only outer coefficients are learned, so training reduces to convex or overdetermined linear least squares (He et al., 30 Mar 2026, Cheng et al., 18 Mar 2026). In the Decoupled-DFNN method for Navier–Stokes, nonlinearity is handled by Gauss–Newton linearization of the stream-function or vector-potential subproblem, after which each iteration again reduces to a linear least-squares solve (Cheng et al., 18 Mar 2026). By contrast, DAF-FlowNet is an unsupervised coordinate-based implicit network trained per timeframe with Adam and a cosine data-consistency loss, while the PeFNO stress models use FNO backbones whose outputs are passed through fixed Fourier-space curl layers (Bisbal et al., 31 Mar 2026, Khorrami et al., 2024).

From a function-space viewpoint, universality is established in several settings. The differential-forms constructions of Neural Conservation Laws are universal for divergence-free vector fields on periodic domains up to harmonic components (Richter-Powell et al., 2022). RTNNs prove a corresponding universal approximation theorem for sufficiently smooth divergence-free symmetric tensors (Jnini et al., 2 Mar 2025). The operator-theoretic viewpoint from orthogonal recurrent models adds a complementary statement: divergence-free dynamics correspond to skew-symmetric directional derivative operators, while soft divergence penalties enlarge the model class to general vector fields (Azencot et al., 2021).

4. Applications and empirical behavior

In incompressible flow solvers, dfNNs are used both as ansatz spaces and as learned dynamics models. The Decoupled-DFNN method for Stokes and Navier–Stokes represents velocity by a stream function in $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$3D or a vector potential in $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$4D, achieves divergence errors at machine precision $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$5–$J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$6 in experiments, and reports approximately $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$7 speedups relative to coupled TransNet solves in several regimes (Cheng et al., 18 Mar 2026). In neural-operator and generative settings, divergence-free projection markedly improves long-term stability: the project-and-generate framework reports exact incompressibility up to discretization error and substantially improved stability and physical consistency on $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$8D Navier–Stokes, while divergence-free diffusion models for Kolmogorov flow show that manifold projection gives $J(t,x) := (\rho(t,x),\rho(t,x)v(t,x)),$9 and that predictor–corrector Leray variants preserve exact incompressibility with competitive rollout errors (Li et al., 25 Mar 2026, Genuist et al., 27 Jan 2026). At the discrete LES level, a divergence-consistent formulation coupled with a CNN closure is stable and accurate with both a-priori and a-posteriori training, whereas divergence-inconsistent filtering is unstable unless stabilized by additional measures (Agdestein et al., 2024).

In solid mechanics, dfNNs encode quasi-static equilibrium directly in the architecture. The physics-encoded FNO (PeFNO) outputs a stress potential rather than stress, and a non-trainable Fourier-space curl mapping enforces $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$0 for all parameters (Khorrami et al., 2024). In reported comparisons, the weak-divergence baseline exhibits pixel-wise errors around $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$1 for $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$2 and $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$3, whereas the divergence-free model yields significantly smaller equilibrium residuals across $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$4 (Khorrami et al., 2024). A later stress-potential formulation reports that the largest error in $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$5 from PeFNO is three orders of magnitude smaller than PiFNO with $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$6, at comparable stress accuracy (Khorrami et al., 1 May 2026). RTNNs broaden this pattern from vector fields to conservative tensor surrogates: for the Cylinder benchmark at $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$7, RTNN reports relative $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$8 error $$\operatorname{div}_{t,x} J = \partial_t \rho + \nabla \cdot (\rho v).$$9 versus $J$0 for Neural Conservation Laws and $J$1 for PINN; for the $J$2D Beltrami case at $J$3, RTNN reports $J$4 versus $J$5 and $J$6 (Jnini et al., 2 Mar 2025).

Medical imaging and geophysics provide application domains in which divergence-freeness encodes mass conservation rather than fluid-state admissibility alone. DAF-FlowNet parameterizes velocity as the curl of a vector potential and combines this with a cosine data-consistency loss that exploits phase periodicity. On synthetic aortic $J$7D Flow MRI, it reports up to $J$8 lower velocity normalized root mean square error, $J$9 lower directional error, and $J^\flat$0 lower divergence relative to the best-performing alternative across noise levels; at peak velocity/velocity-encoding ratios of $J^\flat$1 and $J^\flat$2, it leaves $J^\flat$3 and $J^\flat$4 residual wrapped voxels, respectively (Bisbal et al., 31 Mar 2026). On $J^\flat$5 hypertrophic cardiomyopathy datasets, aortic AAoP–AAoD flow bias is reduced from $J^\flat$6 to $J^\flat$7 in HNCM, and pulmonary bias from $J^\flat$8 to $J^\flat$9 (Bisbal et al., 31 Mar 2026). For Antarctic ice flux interpolation over Byrd Glacier, the stream-function dfNN achieves RMSE $1$0, MAE $1$1, and MAD $1$2; with directional guidance the RMSE improves to $1$3, while PINN retains nonzero MAD and unconstrained NN exhibits large divergence artifacts (Bente et al., 7 Oct 2025).

5. Hard constraints, soft penalties, and common points of confusion

A recurring comparison is between physics-encoded dfNNs and penalty-based models. In the cited works, the hard-constraint approach removes the need to balance a divergence penalty against data fit or PDE residuals. This is emphasized in the finite-neuron approximation theory, in solid-mechanics PeFNO models, in DAF-FlowNet, and in ice-flux interpolation, where soft penalties are described as requiring delicate weighting and as being prone to optimization or stability difficulties (He et al., 30 Mar 2026, Khorrami et al., 2024, Bisbal et al., 31 Mar 2026, Bente et al., 7 Oct 2025). The DNN-MG study provides an explicit negative example: adding penalties on the divergence of the correction or on the corrected field did not significantly reduce divergence, and for $1$4 the network output collapsed toward zero, whereas the stream-function architecture produced better fidelity (Margenberg et al., 2020).

A second point is that exact divergence-freeness is not equivalent to full physical correctness. Curl-based stress parameterizations enforce equilibrium but do not enforce symmetry of the first Piola–Kirchhoff stress, which is not generally symmetric in nonlinear elasticity (Khorrami et al., 2024). In the diffusion setting, exact projection can improve short-horizon fidelity yet be less favorable than softer corrections in some out-of-distribution regimes: for the perturbed Taylor–Green vortex, the PINN-like soft-constraint variant reports the best MSE, while hard-projected variants preserve $1$5 (Genuist et al., 27 Jan 2026). Likewise, in DAF-FlowNet, the cosine loss becomes difficult in the very low VENC regime, where specialized unwrapping methods remain preferable (Bisbal et al., 31 Mar 2026).

A third point concerns representation versus geometry. The theorem linking divergence-free vector fields to skew-symmetric directional derivative operators does not imply that every divergence-free field has the linear form $1$6; the finite-dimensional skew-symmetric matrix is only the linear special case, whereas the general statement is operator-theoretic (Azencot et al., 2021). More generally, coexact or curl representations may require harmonic components on periodic or topologically nontrivial domains, and gauge choices matter whenever the potential is not uniquely determined (Richter-Powell et al., 2022).

6. Domain assumptions, limitations, and extensions

Most exact dfNN constructions rely on structural assumptions about the domain or discretization. The divergence-free finite-neuron theory assumes a bounded Lipschitz domain that is starlike with respect to a ball; on general domains, a finite-dimensional harmonic subspace must be added through Hodge decomposition (He et al., 30 Mar 2026). Spectral Leray projection and the diffusion-based models are presented on periodic domains, with the zero mode removed and FFT-based implementations throughout (Li et al., 25 Mar 2026, Genuist et al., 27 Jan 2026). Stress-potential PeFNO models likewise rely on periodic Fourier representations of unit cells (Khorrami et al., 2024, Khorrami et al., 1 May 2026).

Regularity and conditioning are equally central. In the ReLU$1$7 finite-neuron theory, the approximation rate depends on $1$8, and $1$9 rates require $\mathrm{div}\,P = 0$00; low regularity, such as corner singularities in lid-driven cavity flow, degrades observed convergence (He et al., 30 Mar 2026). The Decoupled-DFNN potential formulation involves fourth-order operators such as $\mathrm{div}\,P = 0$01 or $\mathrm{div}\,P = 0$02, which are handled by automatic differentiation and least squares but can challenge stability at high Reynolds number (Cheng et al., 18 Mar 2026). RTNNs require second derivatives of scalar coefficient fields, and further derivatives appear when viscous or MHD residuals are added (Jnini et al., 2 Mar 2025). DAF-FlowNet, although robust to moderate segmentation perturbations, incurs substantial runtime—$\mathrm{div}\,P = 0$03 on synthetic data and $\mathrm{div}\,P = 0$04 in vivo on an NVIDIA A100, compared with much faster baselines (Bisbal et al., 31 Mar 2026).

The reported extensions are broad but structured. The divergence-free finite-neuron framework explicitly suggests deep networks, adaptive parameter distributions, and analogous treatments of curl-free fields, Maxwell complexes, and neural-operator trial spaces (He et al., 30 Mar 2026). Stress-potential work proposes incorporating nonzero body forces, anisotropy, plasticity, and alternative neural operators such as DeepONets (Khorrami et al., 1 May 2026). DAF-FlowNet points to multi-VENC acquisitions, adaptive Fourier-feature scaling, and broader anatomical validation (Bisbal et al., 31 Mar 2026). Projection-based operator learning identifies non-periodic boundary conditions, $\mathrm{div}\,P = 0$05D vector-potential formulations, and more general geometries as natural next steps (Li et al., 25 Mar 2026, Genuist et al., 27 Jan 2026). Taken together, these directions suggest that dfNNs are less a single architecture than a design principle: encode the relevant divergence constraint in the representation, projection, or operator pipeline, and then optimize only within that constrained function class.