Deep Learning PDE Models (DPMs)

• Deep Learning PDE Models (DPMs) are a diverse family of neural methods that solve, discover, and augment partial differential equations using residual minimization, operator-learning, and hybrid strategies.
• They leverage architectures such as PINNs, Deep Ritz, DeepONet, FNO, and neural ODEs to map solution operators, identify governing equations, and correct incomplete physics.
• These models enable efficient inverse problem solving, uncertainty quantification, and adaptive treatments of high-dimensional or noisy data in scientific computing.

Searching arXiv for the cited DPM literature and related survey context. I’ll look up the specific arXiv papers mentioned for confirmation and context. Deep Learning PDE Models (DPMs) denote a broad class of neural approaches for partial differential equations in which deep networks are used to solve PDEs, learn solution operators, discover governing equations from data, augment incomplete physics, or accelerate related inverse problems. In the surveyed literature, this space includes Physics-Informed Neural Networks (PINNs), the Deep Ritz method, operator-learning architectures such as DeepONet and Fourier Neural Operators, recurrent time-stepping models, sparse-regression discovery pipelines, probabilistic generative surrogates, and hybrid models that couple neural networks to classical or black-box solvers (Tanyu et al., 2022, Freund et al., 2019). This suggests that DPMs are better understood as a family of PDE-constrained learning strategies than as a single architecture.

1. Scope and classification

A standard classification distinguishes between residual-based solution models, variational formulations, operator-learning methods, and inverse or parametric extensions. In the review literature, PINNs approximate $u_\theta(x,t)\approx u(x,t)$ and minimize a composite loss

$$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$

while the Deep Ritz method minimizes an energy functional associated with an elliptic problem; operator-learning methods such as DeepONet and FNO instead learn mappings between function spaces, and inverse formulations either augment the input with parameters $\mu$ or jointly learn $u$ and unknown coefficients (Tanyu et al., 2022).

A distinct usage of the term appears in PDE augmentation. There, the known part of the dynamics is retained and a neural correction is embedded directly in the PDE:

$$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$

with $h_\theta$ learning missing or mis-represented physics such as sub-grid-scale stress or discretization error (Freund et al., 2019).

Family	Canonical role	Representative papers
Residual or variational solvers	Enforce PDE residuals or energies	(Tanyu et al., 2022)
Operator and surrogate models	Learn parametric or functional solution maps	(Tanyu et al., 2022, Wang et al., 2023, Deveney et al., 2019, Longhi et al., 12 Feb 2025)
Discovery and augmentation models	Identify PDE structure or embed neural closures	(Xu et al., 2019, Stephany et al., 2022, Freund et al., 2019)

The same review also treats parametric PDEs and parameter identification problems as first-class DPM use cases, rather than as secondary applications. That is important because many later architectures are explicitly built around $(x,t,\mu)\mapsto u$ or around operator maps from coefficients, initial data, or boundary data to full solution fields (Tanyu et al., 2022).

2. Learned time-stepping and latent-evolution solvers

One prominent line of DPM research makes the analogy to classical numerical time integration explicit. Neural-PDE targets time-dependent PDEs of the form

$$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$

with experiments on the 1D wave equation, 1D/2D heat equation, Burgers’ equation, and a coupled Cahn–Hilliard–Navier–Stokes system. Its core design treats the spatial grid as a sequence and applies a bidirectional LSTM encoder followed by a bidirectional LSTM decoder with dense readout. In place of an explicit finite-difference stencil, it learns

$$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$

so that the LSTM state implicitly carries neighbor information and the weights emulate difference operators and nonlinear fluxes on the fly. The reported training setup uses 30 time steps observed and the next 10 predicted, a 2-layer Bi-LSTM encoder and decoder with 20 hidden units, batch size approximately 25 spatial nodes, and 20 epochs sufficing for convergence. Test $L^2$ errors range from $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$0 down to $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$1, with parameter efficiency of roughly $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$2 parameters versus $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$3–$$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$4 in CNN/MLP baselines to achieve comparable accuracy (Hu et al., 2020).

Finite Difference Nets pursue a related goal from a different angle. They interpret each update $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$5 as a neural analogue of a finite-difference step and train the recurrent block by minimizing either a differential-form residual or an integral-form residual. The framework is tested on a one-way wave equation and on 1D Euler shock-tube problems, uses Conv1D or fully connected layers, and does not require large pre-computed datasets. For the one-way wave test, a 3-layer fully connected network with 100 hidden units per layer and ReLU activations is used; for Euler, one configuration uses a single hidden fully connected layer of 330 units, and another uses three separate fully connected modules of 1024 units each. The reported observations are that integral-form training yields roughly twice the accuracy of differential-form residual training, the training loss decays exponentially over about 300 iterations per time step, shocks and rarefactions are sharply captured with relative $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$6 error below $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$7, and stability is maintained for $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$8 (Chang et al., 2021).

A more recent solver family combines dimensionality reduction with Neural ODEs for parametrized and time-dependent PDEs. The construction introduces an encoder $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$9, a decoder $\mu$0, and a latent ODE

$\mu$1

with autoregressive rollout in latent space. Training combines an autoencoder reconstruction loss, a teacher-forcing latent-dynamics loss, an autoregressive rollout loss, an optional time-generalization loss, and an $\mu$2 latent regularizer. On 1D Advection, 1D Burgers, 2D Shallow-Water, and 2D Molenkamp, the reported outcome is that the DR+Neural ODE model consistently outperforms FNO, cFNO, U-Net, MP-PDE, and Transformer-based operators by $\mu$3–$\mu$4 in normalized RMSE, uses approximately $\mu$5 parameters versus $\mu$6–$\mu$7 in baselines, and achieves inference speed of $\mu$8–$\mu$9 per $u$0 trajectories versus $u$1–$u$2 in competing methods (Longhi et al., 12 Feb 2025).

Across these examples, the common organizing principle is parameter sharing across steps or trajectories. In Neural-PDE this sharing is explicitly compared to classical time-stepping schemes, and in latent Neural ODE models it appears as repeated integration of a single learned vector field (Hu et al., 2020, Longhi et al., 12 Feb 2025).

3. Discovery of governing equations from noisy and limited data

Another major DPM branch is PDE discovery. DL-PDE is a three-stage pipeline: first train a neural surrogate $u$3 on noisy observations, then generate a much larger regular grid of “meta-data,” and finally evaluate derivatives by automatic differentiation and recover the PDE by sparse regression. The candidate library can take the form

$u$4

and the regression problem is

$u$5

The method is tested on groundwater flow, convection–diffusion, Burgers, and KdV. Representative recoveries include $u$6 from $u$7 groundwater samples, $u$8 for convection–diffusion with $u$9 samples, $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$0 for Burgers with $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$1 samples, and $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$2 for KdV with $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$3 samples. The summary states that coefficient errors stay below a few percent even with $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$4–$$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$5 added noise and far fewer than full grid samples (Xu et al., 2019).

PDE-LEARN pushes discovery further by co-training a Rational Neural Network $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$6 and a sparse trainable coefficient vector $$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$7. Its hidden PDE is written through a residual

$$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$8

and optimization couples a data-fit term, a collocation residual term, and an iteratively reweighted least-squares sparsity penalty. The method allows arbitrary data locations, multiple datasets that share the same hidden PDE, and activation functions of the form

$$\partial_t u(t,x)=f_\nu(u,\nabla_x u,\nabla_x^2u)+h_\theta(u,\nabla_xu,\nabla_x^2u),$$9

Representative results include viscous Burgers with $h_\theta$0 noise, where the identified equation is $h_\theta$1; Allen–Cahn at $h_\theta$2 noise, with $h_\theta$3; and a 2D wave equation rewritten as $h_\theta$4, identified at $h_\theta$5 noise as $h_\theta$6. The reported trend is that coefficient errors are usually below $h_\theta$7 up to moderate noise, while failures with too few points or very high noise often still recover a subset of the true terms (Stephany et al., 2022).

A related framework uses neural function approximation followed by differential-operator learning on derivative features. In the Burgers example, the recovered law is

$h_\theta$8

which lies within $h_\theta$9 of the true coefficients. The same study emphasizes that scaling and shifting the inputs correspond to coordinate transformations in the discovered PDE, and it reports a striking ODE reduction: with no spatial derivative inputs, a deep network discovers an ODE surrogate $(x,t,\mu)\mapsto u$0 that matches the Burgers solution well on a short time interval and also captures the bulk temporal dynamics in Swedish weather-station temperature data (Berg et al., 2018).

A common misconception is that data-driven PDE discovery is merely sparse regression on finite differences. The discovery literature above instead relies on smooth neural surrogates, automatic differentiation, meta-data generation, and joint optimization, specifically to address noisy and limited data (Xu et al., 2019, Stephany et al., 2022).

4. PDE augmentation and hybrid solver coupling

The DPM framework in the narrow sense of PDE augmentation embeds a neural network in a PDE expressing the known physics and optimizes over the functional form of the PDE. Its training objective compares trusted data $(x,t,\mu)\mapsto u$1 against the PDE solution $(x,t,\mu)\mapsto u$2 over scenarios and time points, and adjoint PDEs are used to compute

$(x,t,\mu)\mapsto u$3

To reduce cost, the stochastic adjoint method samples one scenario and one time interval per iteration. In large-eddy simulation, the method augments the filtered Navier–Stokes equations by replacing the unresolved term $(x,t,\mu)\mapsto u$4 and discretization errors by $(x,t,\mu)\mapsto u$5. The reported empirical outcome is that the DPM outperforms constant-coefficient and dynamic Smagorinsky models, and also significantly outperforms a priori trained models. For filter width $(x,t,\mu)\mapsto u$6 on a $(x,t,\mu)\mapsto u$7 LES mesh, the DPM reproduces the decay of resolved kinetic energy and the LES energy spectrum with close agreement to filtered DNS; a $(x,t,\mu)\mapsto u$8 Smagorinsky simulation is required to approach the same accuracy. Runtime figures on a single Cray XK7 node are approximately $(x,t,\mu)\mapsto u$9 for $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$0 Smagorinsky, $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$1 for $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$2 DPM with divergence-free projection, and $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$3 for $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$4 DPM without projection (Freund et al., 2019).

A complementary hybrid line addresses settings where the PDE solver is a black box and does not support automatic differentiation. In that setting, a coarse-mesh PDE solver is coupled to a graph neural network that upsamples and corrects the solver output to a fine mesh solution. The trainable quantities are both the GNN weights $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$5 and the mesh coordinates $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$6, with zeroth-order gradient estimators used for $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$7. Three estimators are described: coordinate finite differences, Gaussian random-direction estimation, and a coordinate-Gaussian hybrid. On the reported example, exact differentiation through the solver remains best, but zeroth-order methods outperform a frozen-mesh baseline. The quantitative summary gives test MSE $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$8 for the fully differentiable baseline, $$u_t(x,t)=F\bigl(x,u,\nabla_xu,\nabla_x^2u,\dots\bigr),$$9 for frozen mesh, $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$0 for Coordinate-ZO with $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$1, $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$2 for Gaussian-ZO with $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$3, and $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$4 for Gauss-Coord-ZO with $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$5. Warm-starting the neural network before unfreezing the mesh improves the zeroth-order methods from a plateau around MSE $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$6 to about $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$7 (Ma et al., 2024).

These two directions share a central idea: the neural component is not external post-processing but part of the forward model being optimized. In one case the network is an embedded closure; in the other it is an error-correcting extension around a non-differentiable solver (Freund et al., 2019, Ma et al., 2024).

5. Probabilistic, high-dimensional, and geometry-aware extensions

Probabilistic DPMs treat the solution operator itself as a conditional distribution. PDNO modifies denoising diffusion probabilistic models to learn $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$8, where $$v_i^{n+1}=f_\theta\!\bigl(\hat u^n_{i-r},\dots,\hat u^n_i,\dots,\hat u^n_{i+r}\bigr),$$9 denotes PDE inputs such as coefficients or initial or boundary data. The denoising network is a U-Net with residual blocks and optional multi-head self-attention in the bottleneck, the PDE parameter field is concatenated channel-wise with the noisy state, and inference proceeds by reverse-time sampling from $L^2$0. By construction, the model yields both a point predictor and uncertainty quantification through Monte Carlo samples, summarized by mean relative $L^2$1 error and mean standard deviation. On a 2D elliptic benchmark, PDNO with $L^2$2 million parameters reports MRLE values of approximately $L^2$3 at mesh sizes $L^2$4, compared with FNO at approximately $L^2$5; on a highly non-smooth 1D advection test, PDNO reports approximately $L^2$6 versus FNO at $L^2$7 (Wang et al., 2023).

DeepMartNet addresses high-dimensional Dirichlet boundary value problems and eigenvalue problems for elliptic PDEs through Varadhan’s martingale formulation. If $L^2$8 is the diffusion associated with the elliptic operator $L^2$9, then

$$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$00

is a martingale, and the learning objective penalizes discrete martingale increments along Euler–Maruyama paths, together with boundary and normalization terms. The method is tested on linear and nonlinear Poisson–Boltzmann equations, a Laplace eigenproblem, and a Fokker–Planck eigenproblem. Reported results include relative solution errors of $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$01 for a linear PBE in dimension $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$02, $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$03 for a linear PBE in dimension $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$04, and eigenvalue errors of $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$05 for a Laplace eigenproblem in dimension $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$06 and $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$07 for a Fokker–Planck eigenproblem in dimension $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$08, with runtimes from minutes to tens of minutes on an A100 GPU node (Cai et al., 2023).

Geometry-aware DPMs extend the setting beyond structured grids. One framework solves elliptic PDEs on unknown manifolds represented only by point clouds. Diffusion Maps yields a graph-Laplacian approximation of the manifold operator, and a neural network is trained by minimizing a least-squares loss built from that discrete operator together with any boundary penalties. The paper proves consistency when the network class has infinite width or depth, and for sufficiently wide two-layer networks proves global convergence of gradient descent to a global minimizer of the empirical loss. Numerical examples range from spheres and tori to the Stanford Bunny and a Human-Face geometry; across these tests, the neural solution generalizes to new points with test errors almost identical to training errors and supersedes Nyström interpolation in difficult scenarios (Liang et al., 2021).

Surrogate-based Bayesian inversion provides another extension of DPMs from forward solution to uncertainty-aware parameter inference. A global neural approximation $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$09 is trained over both spatial coordinates and parameters, after which Markov chain Monte Carlo replaces repeated numerical solves by cheap network evaluations. In voltammetry and a heat-transfer problem with a function-valued Biot number, the reported acceleration is $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$10–$$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$11, with surrogate training costs of about $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$12–$$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$13 minutes amortized over large MCMC chains (Deveney et al., 2019).

6. Mathematical structure, limitations, and open questions

The literature attributes several recurring strengths to DPMs. They can exploit automatic differentiation for $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$14, $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$15, and higher derivatives; they can treat parametric PDEs by augmenting inputs with $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$16; operator methods offer fast test-time evaluation after expensive offline training; and augmentation methods can optimize over the full PDE rather than fit local closures in isolation (Tanyu et al., 2022, Freund et al., 2019). Some models are explicitly designed for noisy or sparse measurements, such as DL-PDE and PDE-LEARN, while others focus on uncertainty quantification, such as PDNO, or on dimension-independent stochastic formulations, such as DeepMartNet (Xu et al., 2019, Stephany et al., 2022, Wang et al., 2023, Cai et al., 2023).

The limitations are equally consistent across papers. Neural-PDE notes that no rigorous stability or convergence theory is yet available and that long-term rollout may accumulate error; FD-Nets likewise state that rigorous stability and consistency theory remains to be developed and that on-the-fly optimization at each time step incurs overhead; PDE-LEARN requires a sufficiently rich candidate library, can lose small true coefficients near the pruning threshold, and degrades when $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$17 or noise $$L(\theta)=L_{\text{data}}+L_{\text{PDE}}+L_{\text{bc}},$$18; the black-box hybrid solver finds that zeroth-order differentiation underperforms exact automatic differentiation; and DeepMartNet requires careful tuning of penalty weights and large numbers of SDE paths for high accuracy (Hu et al., 2020, Chang et al., 2021, Stephany et al., 2022, Ma et al., 2024, Cai et al., 2023).

At the theoretical end, a reciprocal line of work formulates deep learning itself through a continuum-limit integro-differential PDE, with forward well-posedness, existence of minimizers for the learning problem, backward co-state equations for stationarity, Pontryagin minimum principles, and an infinite-dimensional Hamilton–Jacobi–Bellman equation for optimal feedback control (Markowich et al., 2024). Although this line is not a PDE solver in the usual sense, it places DPMs within a wider PDE-constrained optimization framework and provides a mathematical language for stability, controllability, and forward-backward discretization.

A plausible implication is that the main distinctions among DPMs are not between “neural” and “numerical,” but between different locations at which PDE structure is imposed: in the loss, in the architecture, in a latent dynamical system, in a stochastic representation, in a sparse symbolic library, or in an embedded correction to a classical solver. That plurality is the defining feature of the DPM literature as represented by current arXiv work (Tanyu et al., 2022, Markowich et al., 2024).