Advancing High-Fidelity Neural Methods: Toward Fast and Accurate Solvers in Scientific Computing and Computational Fluid Dynamics

The accurate numerical solution of partial differential equations is a central problem in scientific computing, at the core of simulation and analysis across engineering, physics, and the applied sciences. Physics-Informed Neural Networks (PINNs) have emerged as a mesh-free alternative to classical discretization methods, yet they fail on the problems that matter most: stiff systems, multi-scale dynamics, and equations with structural conservation laws. This thesis addresses these failures through five contributions, with computational fluid dynamics as the primary application domain and methods evaluated across the full range of PDE problem classes. Training PINNs is notoriously difficult due to the non-convexity of the residual loss landscape; standard first-order optimizers fail to exploit the underlying function-space geometry of the residual operator, leading to ill-conditioned updates and mediocre accuracy even on smooth problems. Gauss–Newton Natural Gradient Descent (GNNG) addresses this by discretising Gauss–Newton’s method in function space, yielding a natural gradient update that provably mimics function-space dynamics and achieves close to single-precision accuracy on the Navier–Stokes benchmarks considered, to our knowledge representing one of the strongest accuracy results reported for PINN-based solvers on these problems. The canonical Gauss–Newton step incurs O(n³) cost in the number of network parameters, making it impractical at scale. By interpreting the Gauss–Newton update through a primal-dual lens, Dual Natural Gradient Descent (D-NGD) derives an equivalent step whose cost scales with the residual space dimension rather than the parameter count; augmented with geodesic acceleration and Nyström preconditioning, this enables second-order PINN training of networks with up to 12.8 million parameters on a single GPU. A systematic benchmark across elliptic, parabolic, hyperbolic, and stiff ODE problem classes establishes that the natural gradient family achieves the highest accuracy on all smooth problems, while confirming that for shock-dominated hyperbolic systems the loss formulation is the binding constraint. Even with an optimal optimizer, standard PINN training demands repeated backpropagation through differential operators at every iteration, dominating wall-clock time in high-accuracy regimes. Gradient-Free Homotopy Residual Learning (GF-HRL) eliminates this cost entirely by freezing the nonlinear feature map and reducing training to a sequence of linear coefficient updates driven by a single precomputed factorisation of a frozen linear operator, cutting training time by two to four orders of magnitude while matching the accuracy of state-of-the-art gradient-based methods. Finally, even a well-trained PINN enforces conservation laws only as soft penalties, with no guarantee of constraint satisfaction. Divergence-free symmetric tensors are a mathematical structure arising naturally in continuum mechanics, unifying mass and momentum conservation under a single tensorial condition; yet no prior neural architecture has exploited this structure as an inductive bias. Riemann Tensor Neural Networks (RTNNs) are the first such architecture, embedding the divergence-free symmetric tensor condition directly into the network by construction and satisfying it to machine precision for all parameter values, achieving the lowest reported errors across Euler, Navier–Stokes, and magnetohydrodynamics benchmarks.