A Systematic Analysis of Automatic Differentiation versus Discretization-based Constraints for Physics-Informed PDE Solvers

• Existing evaluations leave practical guidance underexplored: accuracy alone is insufficient, and efficiency, robustness, stability, and applicability also matter in engineering.Source: §1 · p. 2
• The study compares AD and discretization-based formulations across PDEs with increasing nonlinearity and boundary complexity, without aiming to declare one universally preferable.Source: §1 · p. 3
• Discretization-based methods’ accuracy advantage becomes more pronounced as nonlinearity intensifies, while AD approximation and optimization errors escalate faster in strongly nonlinear regimes.Source: §1 · p. 3
• Its error analysis separates approximation, discretization, and optimization errors, using numerical solutions under identical discretization settings as truncation-error references.Source: §1 · p. 3
• The evaluation spans accuracy, efficiency, and robustness across paradigms and architectures to derive practical solver-selection guidelines and identify training bottlenecks.Source: §1 · p. 3

2 Method

• The framework compares MLP and GNN architectures, AD and discretization-based constraints, and problems with varying nonlinearity and boundary complexity.Source: §2 · p. 3, Figure 1 · p. 4

2.1 Network Architectures

• MLPs process points independently, whereas GNN message passing aggregates neighbourhood information; the study uses GNNs to compare that added capability under controlled conditions.Source: §2.1 · p. 4

2.2 Continuous and discrete constraints

• Discretized constraints impose stable, deterministic truncation errors on predicted solutions, while AD differentiates through the network and couples physical constraints to the training process.Source: §2.2 · p. 5

2.3 Decomposition of errors

• The authors use high-fidelity DNS references and numerical solutions at matching discretization settings, while noting interpolation error from projecting references onto evaluation points.Source: §2.3 · p. 5
• The analysis defines approximation error as architecture bias, optimization error as constraint-embedding and boundary-enforcement bias, and truncation error as an auxiliary reference anchor.Source: §2.3 · p. 5

3.1 Implementation details

• The experiments compare continuous MLP PINNc, discrete MLP PINNd, and discrete GNN PIGNN using general-purpose settings, equal PDE/boundary weighting, and no case-specific priors.Each experiment is repeated five times with different initializations and averaged; near-wall clustering is allowed as a generic practice, but exact shock positions are excluded.Source: §3.1 · p. 7

3.2 Linear to nonlinear Poisson

• As nonlinearity and mesh refinement increase, neural-network errors rise and high-order schemes lose benefits; discrete constraints regularize optimization but help chiefly at lower resolutions.The authors hypothesize that beyond a critical resolution, refinement can hurt more than help, and that this threshold decreases as nonlinearity strengthens.Source: §3.2 · p. 8, Table 2 · p. 8, Table 3 · p. 9
• For simple Poisson cases, MLPs and GNNs have comparable accuracy, but MLPs train more efficiently and stably; GNN loss can more reliably flag poor performance.Residual loss and L2 error may disagree: higher-order schemes can achieve better loss while shifting larger errors into high-frequency regions.Source: §3.2 · pp. 8–10, Figure 4 · p. 10

3.3 Lid-driven cavity Flow

• At Re = 400, PINNc errors rise sharply and then its predictions fail at higher Reynolds numbers, while collocated-grid PINNd and PIGNN errors grow steadily and nearly linearly.PINNc loss can converge normally despite this hidden failure; discrete constraints with global hard constraints reduce optimization difficulties and compensate for truncation error.Source: §3.3 · pp. 11–12, Figure 6 · p. 12

3.4 Backward-facing step Flow

• At 100×50 resolution, PIGNN predicts velocity better than same-resolution DNS, and it captures step-corner flow details more reliably than PINNc or PINNd under the training budget.The authors attribute PIGNN’s robustness to edge-topology information; interpolation artifacts near the corner singularity can obscure global convergence trends across resolutions.Source: §3.4 · p. 13, Table 6 · p. 13, Figure 7 · p. 14

3.5 Hypersonic inviscid cylinder flow

• For the hypersonic cylinder at Ma = 8, extrapolation-based boundary constraints accelerate convergence, unlike their slight error costs in the moderately nonlinear incompressible-flow cases.The network assessment targets local observables because training does not resolve high-fidelity shock positions; the DNS reference uses a coarse-level solution for comparison.Source: §3.1 · p. 6, §3.5 · p. 15
• PIGNN recovers major expansion features and reasonable wall-surface trends in the hypersonic cylinder, while PINNc fails to propagate shocks and PINNd produces smeared, physically deficient profiles.PINNc is excluded from local-profile comparisons because its results are insufficient quality; PIGNN remains indicative for preliminary assessment, not high-fidelity engineering use.Source: §3.5 · pp. 15–16
• MUSCL reconstruction and other low-dissipation fluxes destabilize training in the strongly nonlinear hypersonic case, whereas Rusanov flux is more robust and AD-based solvers are still more unstable.MUSCL models remain sub-optimal after 30000 epochs, indicating instability associated with high-order-like neural training in stiff problems.Source: §3.5 · pp. 15–16
• At 50000 epochs, PINNd’s lower boundary and PDE residuals accompany severely over-smoothed shocks, while prolonged training can degrade its predictions without an obvious loss warning.This illusory convergence also occurs under discretization constraints and may reflect intrinsic MLP optimization; PIGNN stays near its achievable bound with stable loss oscillations.Source: §3.5 · pp. 16–17
• PIGNN remains stable through Mach 25, with predictions close to theoretical values, but viscous-flow attempts fail amid near-wall and shock-layer smearing.Higher Mach cases require 20000–40000 epochs; at elevated Mach, the temperature peak shifts slightly upstream while shock location does not appreciably move forward.Source: §3.5 · p. 17

3.6 Efficiency analysis

• Neural-solver efficiency must include preprocessing and robustness: PINNc’s mesh-free advantage applies to weakly nonlinear cases, while mesh-based solvers offer simpler collocated-grid workflows than staggered-grid implementations.DNS remains important for final high-fidelity evaluation because its error sources are traceable and controllable; extra neural-solver cost with limited performance gain is undesirable.Source: §3.6 · pp. 18–19
• PIGNN’s runtime grows faster than PINNc’s with mesh refinement for the cavity and backward-facing step, while in hypersonic flow it trains stably at a time comparable to explicit schemes.The hypersonic case has 6480 cells and four outputs; FV residual assembly uses face fluxes and one sparse scatter-add, while FD constraints require multiple matrix operators and block arithmetic.Source: §3.6 · pp. 18–19, Table 7 · p. 19

4 Conclusion and discussion

• AD performs well in weakly nonlinear, well-conditioned regimes with globally enforceable hard boundary constraints, but becomes increasingly vulnerable to illusory convergence and abrupt degradation as nonlinearity strengthens.Nonlinearity strength, including source-term stiffness or Reynolds number, is proposed as a transferable criterion for locating this boundary.Source: §4 · pp. 19–20
• For discretization-based formulations, stronger nonlinearity shifts dominant error from truncation toward approximation and optimization, while higher-order schemes can lose or reverse their benefits.MSE averaging also limits attainable accuracy: mesh-refinement gains can stagnate or degrade beyond a threshold that decreases with nonlinearity.Source: §4 · p. 20
• GNNs remain more robust under strong nonlinearity and complex boundary conditions, whereas their aggregation capacity adds computational overhead without appreciable accuracy gains on simple weakly nonlinear problems.Their loss trajectories are more diagnostically reliable than MLP trajectories, which can show deceptive convergence without visible precursors.Source: §4 · p. 20
• Extrapolation-type boundary conditions act like deterministic truncation-error perturbations, adding marginal error in weakly nonlinear problems but accelerating convergence in strongly nonlinear, ill-conditioned cases.The paper therefore treats them as problem-dependent regularization rather than universally beneficial or harmful.Source: §4 · p. 20
• Engineering evaluation should delineate each method’s applicability and performance ceiling, and weigh marginal accuracy gains against computational costs.The paper also calls for formulation-agnostic diagnostics that remain informative when loss histories are misleading, reducing reliance on unavailable a-posteriori reference solutions.Source: §4 · p. 20

A Details of numerical solvers

• All numerical discretizations avoid conventional ghost points and ghost cells to remain consistent with the neural solvers’ constraint formulations.The illustrated layouts are Cartesian-grid standard FD, polar-coordinate FD, staggered-grid FD, and vertex-centered FV.Source: A Details of numerical solvers · p. 23
• Figure A1 depicts Cartesian nodal, polar-coordinate nodal, staggered-grid variable, and vertex-centered FV dual-control-volume arrangements.The detailed discretization strategies are described in the following subsections.Source: A Details of numerical solvers · p. 23
• The appendix presents these layouts as the numerical arrangements used for the corresponding discretizations.The schematic labels Cartesian and polar nodal grids, staggered-grid variables, and a dual control volume for vertex-centered FV.Source: A Details of numerical solvers · p. 23

A.1 FDM for Poisson-type equations

• The fourth-order Cartesian Laplacian stencil uses two neighboring nodes on each side and has truncation error O(∆x^4) + O(∆y^4).Near boundaries, the fourth-order scheme reduces to second-order stencils; Dirichlet conditions are imposed at boundary nodes.Source: §A.1 · pp. 23–24
• Cartesian Poisson discretizations produce a sparse system Au = f, solved directly for the linear Sine case and by Newton iteration for the nonlinear Polynomial case.The appendix presents second- and fourth-order central-difference approximations for the Laplacian.Source: §A.1 · pp. 23–24
• Polar-coordinate Liouville discretization uses second- and fourth-order central differences, with second-order stencils near the outer boundary and polar center.Because direct polar stencils are ill-defined as r → 0, the polar-center value is discretized by averaging neighboring nodes on the first ring.Source: §A.1 · p. 24
• The stiff Liouville source-term problem uses higher-order Halley iteration because conventional Newton iteration is unstable.The Halley update approximates the quadratic-term increment using a Newton step; practical iterations still require relaxation factors.Source: §A.1 · pp. 24–25

A.2 Staggered-grid FDM for incompressible NS equations

• The staggered grid stores horizontal and vertical velocity at their respective half-grid nodes and pressure at primary nodes.The arrangement addresses pressure–velocity decoupling associated with conventional collocated grids.Source: §A.2 · p. 25
• At moderate Reynolds numbers, the cases use standard central differences and half-grid difference operators for pressure gradients and continuity-equation velocity gradients.The half-grid operators are equivalent to standard second-order finite-difference approximations.Source: §A.2 · p. 25
• The Navier–Stokes equations are cast as a nonlinear algebraic system G(U) = 0, with separate discrete operators for velocity gradients, pressure gradients, viscous terms, and pressure constraints.The solution vector is U = [u, v, p]ᵀ, and boundary-condition treatments are embedded in the viscous-term matrices.Source: §A.2 · p. 25
• Because velocity is not stored at physical boundaries on the staggered mesh, quadratic-fitting polynomial constraints impose velocity boundary conditions.The prescribed physical boundary conditions also leave the implicit discrete system rank-deficient, requiring supplementary Neumann-type conditions.Source: §A.2 · p. 25
• Newton iteration with an analytical Jacobian solves the nonlinear algebraic system, and GMRES solves each linear subproblem.The staggered-mesh solutions are mapped to vertex nodes by linear interpolation.Source: §A.2 · pp. 25–26

A.3 Vertex-centered FVM for compressible Euler equations

• The steady system is advanced by pseudo-time marching in conservative variables, using first-order forward Euler steps initialized from the free-stream state.The solution vector is U = [ρ, ρu, ρv, ρE]ᵀ, and the spatial operator collects the inviscid flux terms.Source: §A.3 · p. 26
• The formulation uses non-dimensional free-stream reference conditions, with the cylinder radius as characteristic length and all free-stream reference quantities set to unity.Variables are treated as non-dimensional after the scaling relations are applied.Source: §A.3 · p. 26
• The differential equations are converted to integral form over vertex-centered dual cells, and edge fluxes are aggregated into nodal residual contributions.Dual-cell areas are exact on orthogonal meshes and use an area-averaging approximation on non-orthogonal meshes.Source: §A.3 · pp. 26–27
• The method reconstructs left and right dual-edge states with MUSCL and a Venkatakrishnan limiter, then uses Rusanov fluxes to form residuals.The limiter uses the minimum edge-wise coefficient at each vertex; the Rusanov flux combines a central term with maximum-wave-speed dissipation.Source: §A.3 · p. 27
• Boundary constraints are converted from primitive to conservative variables and enforced iteratively, while wall temperature is used to check convergence because residual oscillations persist.The hypersonic case does not achieve the clean convergence behavior seen in incompressible tests; iteration stops when wall-temperature variations cease to be notable.Source: §A.3 · p. 27

B BC enforcement and loss functions

• Training is unsupervised, with a composite loss combining PDE residual and boundary-condition MSE terms; fully hard-constrained cases set the boundary loss to zero.Partially hard-constrained cases use equal weights for the two terms, except hypersonic cases with additional weighting for magnitude discrepancies.Source: B BC enforcement and loss functions · pp. 27–28

B.1 Poisson-type equations

• For Poisson-type cases, global hard boundary constraints use a distance-function ansatz, while discrete models impose algebraic residuals at interior nodes.The network prediction and hard-constraint treatments are shared across baseline models; the PDE residual formulation differs for discretization-based models.Source: §B.1 · p. 28

B.2 Lid-driven cavity flow

• For the lid-driven cavity, collocated-grid discrete operators represent gradients and viscous terms, while equal weighting is used for momentum and divergence residuals.A boundary-zeroing factor in the trial solution enforces global hard boundary constraints, so the boundary loss is identically zero.Source: §B.2 · p. 28

B.3 Backward-facing step flow

• Global hard boundary constraints are unavailable for the step, so constraints are imposed only on partial boundaries with a unified boundary-residual formulation.Terms on hard-constrained segments vanish identically, and the stated ansatz fully satisfies the pressure outlet condition and selected velocity constraints.Source: §B.3 · pp. 28–29
• For discrete models PINNd and PIGNN, DNS-compatible outlet extrapolation constraints may be enforced explicitly or omitted, with neighboring boundary nodes coupled through discrete stencils.The table reports global relative L2 errors across outlet treatments and mesh resolutions at Re = 100.Source: §B.3 · p. 29, Table B1 · p. 29
• Omitting outlet extrapolation constraints yields slightly lower velocity errors but larger pressure errors, while enforcing them reverses that trend; PINNd is more sensitive than PIGNN.The authors omit these constraints because pressure-gradient quantities are of greater practical interest than absolute pressure in incompressible engineering problems.Source: §B.3 · p. 29, Table B1 · p. 29

B.4 Hypersonic inviscid cylinder flow

• Density and pressure are constrained to be positive through a softplus transformation, and the network predicts primitive variables directly for the steady Euler equations.The formulation does not require conservative variables as in pseudo-time-marching solvers.Source: §B.4 · p. 29
• The continuous PINNc uses differential PDE residuals, whereas discrete models use integral Rusanov face fluxes; residual components receive fixed weights based on Mach-number scaling.Momentum and energy fluxes scale as Ma and Ma^2 relative to mass flux, respectively, motivating the component weighting.Source: §B.4 · p. 29
• The inlet coordinate is normalized to improve training stability, and hard Dirichlet conditions are imposed using the free-stream physical values.The characteristic length is specified as 3 for the coordinate normalization, with the free-stream state values given in the implementation passage.Source: §B.4 · p. 30
• Outlet points are included in the continuous PDE residual, and soft boundary penalties impose slip at solid walls; additional zero-normal-gradient conditions are also used.The text reports faster convergence when the zero-normal-gradient outlet condition is additionally imposed.Source: §B.4 · p. 30

C Soft BC enforcement results

• Soft BC penalties produced evidently larger errors than pure hard constraints across the Sine, Polynomial, Liouville, and Re = 100 lid-driven cavity tests.The comparisons used identical training configurations.Source: C Soft BC enforcement results · p. 30, Table C2 · p. 30, Table C3 · p. 31
• PIGNN errors were substantially larger than MLP-based PINNc and PINNd errors, but the gap narrowed from over one order of magnitude on Sine to stronger coupled Navier–Stokes cases.PIGNN also showed higher parameter sensitivity and could diverge at high resolution; architecture changes, hyperparameter tuning, or longer training could mitigate this.Source: C Soft BC enforcement results · p. 30
• For MLP models, AD and high-order discrete operators brought no obvious gain on Sine and Polynomial, while pure hard constraints substantially suppressed interference among optimization objectives.Source: C Soft BC enforcement results · p. 30
• A slash marks configurations with frequent optimization instability and no consistent convergence across repeated trials.Source: C Soft BC enforcement results · p. 30
• The soft-penalty experiments were designed for identical-setting comparison, not to represent the best achievable soft-BC performance.The study also avoided ground-truth-guided post-hoc parameter tuning.Source: C Soft BC enforcement results · pp. 30–31

D Hyperparameter analysis

• For comparable parameter budgets, deeper, narrower PINNc architectures outperformed shallower, wider ones in the hyperparameter analysis.The study tested PINNc configurations once; after network size passed a threshold, further depth or width increases yielded only marginal gains.Source: D Hyperparameter analysis · pp. 31–32, Figure D2 · p. 32
• Higher-resolution settings generally converged more slowly and, except for Sine, did not surpass lower-resolution accuracy even with extended training.For incompressible Navier–Stokes cases, average L2 error continued falling mainly through pressure convergence while velocity error stagnated or rose slightly.Source: D Hyperparameter analysis · p. 32, Figure D3 · p. 32
• The authors attribute high-resolution underperformance to spectral bias in MSE-based optimization rather than insufficient training duration.They advise prioritizing network depth and note that lower-resolution configurations tend to be more robust as nonlinearity increases.Source: D Hyperparameter analysis · p. 32
• The largest training-duration gains occurred within 10000–20000 epochs; further training reduced error only slightly.Each case and resolution was trained up to 50000 epochs with the architecture held fixed.Source: Figure D3 · p. 32, D Hyperparameter analysis · p. 32