Blending Neural Operators and Relaxation Methods in PDE Numerical Solvers

TL;DR

The paper addresses the challenge of efficiently solving PDEs at scale by combining neural operator regression with classical relaxation methods. It introduces HINTS, a hybrid iterative framework that learns a solution operator $\mathcal{G}: (k,f) \mapsto u$ via DeepONet and couples it to standard relaxers by alternating updates in a ratio $1:(n_r-1)$, enabling balanced convergence across eigenmodes and discretization transfer between $\Omega^{h_D}$ and $\Omega^{h}$. The authors demonstrate the approach on Poisson and Helmholtz problems in 1D, 2D, and 3D, including irregular geometries, and extend it with HINTS-MG for multiscale solvers and as a preconditioner in Krylov methods; the framework generalizes across discretizations and domains without retraining. The results show substantial speedups, robustness to indefiniteness, and scalability to large systems, suggesting practical impact for engineering simulations and HPC workflows where traditional solvers struggle or are expensive.

Abstract

Neural networks suffer from spectral bias having difficulty in representing the high frequency components of a function while relaxation methods can resolve high frequencies efficiently but stall at moderate to low frequencies. We exploit the weaknesses of the two approaches by combining them synergistically to develop a fast numerical solver of partial differential equations (PDEs) at scale. Specifically, we propose HINTS, a hybrid, iterative, numerical, and transferable solver by integrating a Deep Operator Network (DeepONet) with standard relaxation methods, leading to parallel efficiency and algorithmic scalability for a wide class of PDEs, not tractable with existing monolithic solvers. HINTS balances the convergence behavior across the spectrum of eigenmodes by utilizing the spectral bias of DeepONet, resulting in a uniform convergence rate and hence exceptional performance of the hybrid solver overall. Moreover, HINTS applies to large-scale, multidimensional systems, it is flexible with regards to discretizations, computational domain, and boundary conditions.

Figures (14)

• Figure 1: Overview of the Hybrid Iterative Numerical Transferable Solver (HINTS). The goal is to solve $u(\mathbf{x})$ in the differential equation. HINTS starts by discretizing the computational domain $\Omega$ into $\Omega^h$, where $h$ denotes the mesh size. We then initialize a guess of the solution ($\mathbf{v}^h=0$, for example). In the $k_\text{it}$th iteration, the approximate solution $\mathbf{v}^h$ is corrected by $\delta\mathbf{v}^h$ through either the DeepONet solver or the standard numerical solver. The choice of the solver is determined by whether $n_\text{r}$ divides $k_\text{it}$, where $n_\text{r}$ is a parameter that we choose to optimize the performance. There is a plurality of choices for the numerical solver, such as Jacobi, Gauss-Seidel, and multigrid methods. The algorithm proceeds until the residual of the underlying PDE is smaller than a threshold $\varepsilon$.
• Figure 2: Results of 1D Poisson Equation. (A) Profiles of $k(x)$ and $f(x)$. (B) Eigenmodes $\bm{\phi}^h_i$ ($\phi_i(x)$; $i=1,5,10$ herein) and corresponding loading vectors $\frac{\mathbf{A}^h\bm{\phi}^h_i}{||\mathbf{A}^h\bm{\phi}^h_i||}$. Mode 1 has the lowest spatial frequency, while mode 10 has a relatively high spatial frequency. (C) Numerical results. We consider three setups, each shown in one row: (M1) Jacobi solver only; (M2) Jacobi solver with DeepONet initializer, i.e., one-time usage of DeepONet followed by Jacobi iterations; (M3) HINTS-Jacobi (with a DeepONet-to-Jacobi ratio $1:24$). The second column shows key snapshots of the iterative solution. The third column shows the histories of the norms of residual and error of the approximate solution, with the snapshots in the second column marked correspondingly. The fourth column shows the history of the norm of error for eigenmodes 1, 5, and 10.
• Figure 3: Results of 1D Helmholtz Equation. (A) Profiles of $k(x)$ and $f(x)$. (B) Eigenmodes $\bm{\phi}^h_i$ ($\phi_i(x)$; $i=1,5,10$ herein) and corresponding loading vectors $\frac{\mathbf{A}^h\bm{\phi}^h_i}{||\mathbf{A}^h\bm{\phi}^h_i||}$. Mode 1 has the lowest spatial frequency, while mode 10 has a relatively high spatial frequency. (C) Numerical results. We consider three setups, each shown in one row: (M1) Jacobi solver only; (M2) Jacobi solver with DeepONet initializer, i.e., one-time usage of DeepONet followed by Jacobi iterations; (M3) HINTS-Jacobi (with a DeepONet-to-Jacobi ratio $1:14$). The second column shows key snapshots of the iterative solution. The third column shows the histories of the norms of residual and error of the approximate solution, with the snapshots in the second column marked correspondingly. The fourth column shows the history of the norm of error in eigenmodes 1, 5, and 10.
• Figure 4: Results of HINTS-Jacobi for 2D Poisson Equation and 3D Helmholtz Equation. The DeepONet-to-Jacobi ratio is $1:14$ in both cases. (A) 2D Poisson equation. In the first row, the first and the second columns shows the discretization of DeepONet and the numerical problem, respectively; the third and the fourth columns show the histories of norms of errors and residuals and mode-wise errors, respectively. In the second row, we display the profiles of errors of some key snapshots together with the true/converged solution. (B) 3D Helmholtz equation. In the first row, the first column illustrates the 3D domain where we define the problem. Other displays are similar to Poisson equation in (A). The snapshots show the errors and true solution on cross section $z=0.5$, marked in the first panel.
• Figure 5: Results of HINTS-MG for 1D Poisson Equation and 2D Helmholtz Equation. (A) 1D Poisson Equation. The first column shows the setup of V-cycle in the multigrid algorithm, where each dot represents 10 iterations. The second and the third column compare the histories of norms of errors and residuals for the two setups of relaxation: (M1) with Gauss-Seidel relaxation only; (M2) with hybrid DeepONet-GS relaxation ($1:9$ in each grid level). (B) 2D Helmholtz equation. Display is similar to (A). The DeepONet-GS ratio is also $1:9$.