1 Introduction↩︎

The Helmholtz equation, which governs time-harmonic wave propagation, plays a fundamental role in a wide range of applications, including acoustics, electromagnetics, geophysics, and medical imaging. Despite its importance, computing accurate numerical solutions remains notoriously challenging, particularly in the presence of heterogeneous media and high wavenumbers. At high wavenumbers, standard discretizations suffer from the pollution effect, with a wavenumber-dependent phase error that accumulates across the domain [1]. In multilevel settings this difficulty is amplified: using the same fixed stencil on progressively coarser grids often yields larger dispersion errors, leading to phase mismatches between levels [2]. Meanwhile, the discretized Helmholtz operator is strongly indefinite and may be close to singular in regimes of practical interest, which breaks the assumptions behind classical multigrid. Standard relaxations then cease to be reliable smoothers, and coarse-grid correction may no longer approximate the dominant error components. As a result, multigrid iterations can stagnate or diverge, with the failure often more pronounced in heterogeneous media [3].

These challenges have motivated a broad class of Helmholtz-specific solvers and preconditioners. Within the multigrid family, shifted-Laplacian methods [4] use a complex-shifted operator to add damping, improving spectral properties and yielding effective Krylov preconditioners. Wave-ray multigrid [5], [6] augments standard multigrid with ray-based corrections to represent oscillatory components that cannot be resolved on very coarse grids. Beyond multigrid, the Convergent Born Series (CBS) [7] constructs a modified Born expansion based on a damped reference medium and a suitable preconditioner, yielding a provably convergent iteration even for large, strongly scattering inhomogeneous media. Domain decomposition methods [8], [9] couple local subdomain solves through inter-subdomain transfer mechanisms, enabling efficient sweeping-type solvers and scalable implementations. Despite these advances, achieving robustness at high wavenumbers and in strongly heterogeneous media without sacrificing efficiency remains an open problem [3].

A useful way to state the multigrid obstruction is that scalar coarsening is not well matched to the phase-space structure of high-frequency waves. After relaxation, a Helmholtz error is not generally smooth in \(x\); in the high-frequency regime it is more naturally described by a local WKB/geometrical-optics ansatz \(a(x)e^{i\phi(x)}\), with wave vector \(\nabla\phi(x)\) near the characteristic shell \(|\xi|\approx k(x)\) [10], [11]. Coarsening only the physical coordinate \(x\) therefore discards local direction, phase, and scattering information that remains essential for correction. This suggests a different coarse object: a coarse node should carry several local wave-packet or phase-space coefficients, not a single scalar degree of freedom [12].

Recently, there has been growing interest in leveraging deep learning for solving partial differential equations (PDEs). Physics-Informed Neural Networks (PINNs) [13] incorporate the governing equations into the training objective and have been applied to a variety of forward and inverse problems, but they typically require retraining for each new parameter configuration and thus generalize poorly across parametrized PDE families. In contrast, neural operator methods [14] aim to learn mappings between infinite-dimensional function spaces. Representative architectures include DeepONet [15], which combines a branch network encoding the input function and a trunk network encoding spatial coordinates, and the Fourier Neural Operator (FNO) [16], which performs global spectral convolutions to capture long-range dependencies. Subsequent extensions further broaden applicability, such as Geo-FNO [17] for irregular geometries, and MgNO [18] with a multigrid-inspired hierarchical design; attention-based operator learning has also been explored, e.g., Galerkin-type attention [19].

Although neural operator methods overcome the generalization limitations of PINNs and enable direct inference, they still face challenges in achieving high accuracy and ensuring solution reliability. This issue stems from the well-known spectral bias of neural networks [20], [21], which tend to learn low-frequency components efficiently while struggling to capture high-frequency modes. To address this, recent studies have proposed Neural Solvers or Preconditioners [22]–[29]. Instead of acting as end-to-end predictors, these methods embed neural networks as learned components within classical iterative frameworks, allowing for progressive error reduction. Notable examples include WAVE-ADR-NS [25], an efficient multigrid algorithm inspired by the Wave-Ray method [6] that learns to eliminate “characteristic” error components specific to the high-wavenumber Helmholtz equation. Similarly, MGCFNN [26] constructs a hierarchical AI solver by combining a neural multigrid architecture with a FNO at the coarse level. A closely related family is hybrid methods [24], [30], [31], which interleave classical iterative solvers with neural operators, using the network to capture the low-frequency error that the conventional smoother struggles to eliminate and thereby accelerate its convergence. While promising, many neural solvers for Helmholtz equations still face practical limitations. Although they reduce the iteration count relative to classical methods, their accuracy and convergence often degrade at high wavenumbers and in strongly heterogeneous media. Moreover, the computational cost of evaluating relatively heavy networks at each iteration can offset these gains, leaving the overall solving time as slow as, or even slower than, that of classical methods.

In this work, we develop Multi-channel Multigrid (McMg), a learned phase-space multigrid preconditioner for heterogeneous Helmholtz equation. The method coarsens physical space as in a standard multigrid hierarchy, but keeps a finite-dimensional representation of unresolved phase, direction, and scattering information in the channel dimension. Our objective is to address a key limitation of existing neural preconditioners: whether efficient solution performance can be retained for high-frequency, high-contrast, and large-scale Helmholtz problems, as required in realistic industrial applications. Our main contributions are summarized as follows:

• We introduce McMg, a learned phase-space multi-channel multigrid preconditioner for heterogeneous Helmholtz problems. The method is formulated as an iterative error-correction scheme rather than an end-to-end neural predictor, and learns an approximation to the inverse operator that can be used as a solver or preconditioner.

• We propose multi-channel multigrid operators based on locally adaptive stencils. The stencil coefficients are generated from the local wave speed, while the residuals and corrections are lifted to a multi-channel space. The channel dimension acts as a learned local packet basis, helping preserve phase and directional information that scalar coarse grids lose.

• We develop a setup strategy with parameter sharing within the multigrid hierarchy. Since the medium coefficient is fixed during the solve, the nonlinear physical features are computed once in the setup phase and cached; every subsequent iteration then applies the cached adaptive operators to the evolving residuals and corrections as a sequence of linear convolutions. Thus, for a fixed medium, the learned V-cycle is linear in the residual.

• We demonstrate cross-scale generalization and large-domain scalability. Models trained on small domains or patches can be deployed on substantially larger domains and higher effective wavenumbers. For further scaling, we introduce Layer-by-Layer Progressive Finetuning (LLPF): by freezing previously trained levels and finetuning only a single newly appended coarse level, LLPF extends the long-range support of the solver to large domains at a fraction of the cost of full retraining, while matching its convergence.

• We validate McMg on high-frequency, high-contrast, and large-scale Helmholtz problems, spanning the OpenBreastUS, OpenFWI CurveFault-B, and Kimberlina 3D benchmarks. Relative to the fastest classical baseline on each benchmark—selected among GMRES with shifted-Laplacian multigrid, CBS, and sparse LU—McMg reduces iteration counts by \(110.6\times\), \(51.3\times\), and \(45.1\times\), with corresponding wall-clock speedups of \(10.3\times\), \(3.4\times\), and \(6.8\times\). Spectral-error, coarse-basis, channel-SVD, and Green’s-function-support diagnostics provide evidence for the learned phase-space interpretation.

2 Preliminaries↩︎

In this section, we establish the mathematical formulation of the heterogeneous Helmholtz equation, discuss the absorbing boundary layers required for domain truncation, and review its finite difference discretization. We then outline the V-cycle Geometric Multigrid algorithm.

2.1 Heterogeneous Helmholtz equation↩︎

We consider the Helmholtz equation, which models wave scattering phenomena in an unbounded heterogeneous medium. To ensure the well-posedness of the solution, we assume that the equation satisfies the Sommerfeld radiation condition. The complete formulation is given by: \[\label{eq:helm} \begin{align} -\Delta u(x) - k(x)^2 u(x) &= f(x), \quad x \in \mathbb{R}^{D},\\ \lim_{|x| \to \infty} |x|^{\frac{D-1}{2}} \left( \frac{\partial u}{\partial |x|} - i k_0 u \right) &= 0, \end{align}\tag{1}\] where \(u(x)\) denotes the wavefield, \(f(x)\) is the source term, \(D\) represents the spatial dimension, and \(i = \sqrt{-1}\) is the imaginary unit. The spatially varying wavenumber is given by \(k(x) = \frac{\omega}{c(x)}\), where \(\omega\) is the angular frequency and \(c(x)\) is the speed of sound. In the radiation condition, \(k_0\) denotes the constant background wavenumber, assuming \(k(x) \to k_0\) as \(|x| \to \infty\).

To facilitate numerical computation, it is necessary to truncate the unbounded domain. However, applying Dirichlet or Neumann boundary conditions on the artificial boundary typically leads to spurious reflections, which do not align with the physical behavior of wave propagation in unbounded domains. To mitigate such reflections, absorbing boundary layers are employed; common techniques include the sponge layer [32] and the perfectly matched layer (PML) [27], [33].

The sponge layer adds a damping term to the original equation, \[-\Delta u(x) - k(x)^2 (1-\gamma i)u(x) = f(x), \quad x \in \Omega \cup \Omega_{\mathrm{abl}},\] where \(\Omega\) denotes the computational domain and \(\Omega_{\mathrm{abl}}\) the absorbing boundary layer; the damping coefficient \(\gamma\) vanishes inside \(\Omega\) and increases polynomially with distance from the interface within \(\Omega_{\mathrm{abl}}\). The PML instead extends the spatial coordinates into the complex plane, introducing artificial decay without reflection through a complex stretching of the Laplacian; we refer the reader to Appendix 6.1.1 for the detailed formulation.

In this work, the computational domain is chosen as a rectangular box (a cube in the three-dimensional case) and discretized on a structured uniform grid of spacing \(h\). We adopt the finite difference method (FDM) as the underlying discretization; although our method ultimately learns adaptive stencils, it builds upon this classical scheme. Due to the heterogeneous wavenumber \(k(x)\), the resulting discrete operator acts as a spatially varying convolution: \[\label{eq:discrete95fdm} \mathcal{K}_{A} * \boldsymbol{u} = \boldsymbol{f},\tag{2}\] where \(\mathcal{K}_{A}\) denotes the local stencil, \(*\) the spatially variant convolution, \(\boldsymbol{u}\) the numerical solution, and \(\boldsymbol{f}\) the discrete source term. In two dimensions, the standard five-point stencil at grid index \((i,j)\) takes the form: \[\label{eq:FDM95stencil} \mathcal{K}_A[i,j] = \frac{1}{h^2} \begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 - (\boldsymbol{k}_{i,j} h)^2 & -1 \\ 0 & -1 & 0 \end{bmatrix},\tag{3}\] where \(\boldsymbol{k}_{i,j}\) denotes the discrete wavenumber \(\boldsymbol{k}\) at location \((i,j)\). Other discretizations, such as the Fourier spectral method (Appendix 6.1.2), are equally applicable. Regardless of the scheme, the resulting linear system can be written in the general form \[A \boldsymbol{u} = \boldsymbol{f}.\]

2.2 Geometric Multigrid (GMG)↩︎

Geometric multigrid (GMG) is an optimal-complexity (\(\mathcal{O}(N)\)) iterative solver for elliptic PDEs that combines fine-grid smoothing with coarse-grid correction over a hierarchy of nested grids [34]. Its efficiency stems from a complementary division of labor: inexpensive local smoothers (such as weighted Jacobi or Gauss–Seidel) rapidly damp the high-frequency components of the error, while the remaining smooth, low-frequency error is transferred to a coarser grid, where it is cheaper to resolve and, relative to the coarser mesh, once again appears oscillatory. Applying this principle recursively yields the V-cycle: starting from the finest grid, the residual is smoothed and then restricted down the hierarchy to the coarsest level \(l=L\), where the system is solved directly; the resulting correction is prolongated back up, with smoothing applied at each level along the way.

While highly effective for elliptic operators, standard GMG can deteriorate or diverge for high-frequency Helmholtz problems [3]. Applying the same fixed stencil across levels yields inconsistent dispersion on coarser grids, so coarse-grid corrections may arrive with an incorrect phase and amplify the error instead of reducing it.

3 Multi-channel Multigrid↩︎

In this section, we detail the proposed Multi-channel Multigrid (McMg) architecture. Unlike classical geometric multigrid methods that operate on scalar fields across levels, McMg uses a multi-channel latent representation. The spatial grid is coarsened, but each coarse node carries several latent coefficients that can encode local phase, propagation direction, and scattering information. Throughout the hierarchy, the discrete operators and smoothers are modulated by the local wave speed, so that the solver adapts to heterogeneous media while preserving wave information often lost in standard scalar GMG.

3.1 Learnable Setup Phase↩︎

The operators used by McMg share a common structure. Given the medium coefficient \(\boldsymbol{m}\) (e.g., the discrete wavenumber \(\boldsymbol{k}\)), a learnable nonlinear map first extracts a physical feature \[\label{eq:phys95feature} \boldsymbol{c} = \sigma(\mathcal{K}_{\boldsymbol{m}} \ast \boldsymbol{m}),\tag{4}\] which then modulates a linear convolution of the field \(\boldsymbol{x}\) (a residual or correction): \[\label{eq:pyhs95conv} \Phi_{\boldsymbol{m}}(\boldsymbol{x}) = \boldsymbol{c} \odot (\mathcal{K}_{\boldsymbol{x}} \ast \boldsymbol{x}).\tag{5}\] Here \(\mathcal{K}_{\boldsymbol{m}}\) and \(\mathcal{K}_{\boldsymbol{x}}\) are learnable multi-channel convolution kernels, \(\ast\) denotes multi-channel convolution, \(\odot\) is element-wise multiplication, and \(\sigma\) is a nonlinear activation (GELU [35] in this work). The construction is deliberately asymmetric in its two inputs. It is nonlinear in the coefficient \(\boldsymbol{m}\), so that the physical feature \(\boldsymbol{c}\) can represent the complex dependence of optimal discretization parameters on the wavenumber (e.g., the dispersion-optimized stencil coefficient \(2\cos(\boldsymbol{k}h)\), which is transcendental in \(\boldsymbol{k}\) rather than a low-order polynomial). For a fixed medium, however, \(\Phi_{\boldsymbol{m}}\) is linear in the field \(\boldsymbol{x}\), matching the linearity of the Helmholtz solution operator \(A^{-1}\): the feature \(\boldsymbol{c}\) acts as a spatially adaptive stencil whose weights are set by the medium and applied linearly to \(\boldsymbol{x}\).

A key consequence of this split is that the physical feature \(\boldsymbol{c}\) in Eq. 4 depends only on the medium coefficient, which is fixed for a given problem. It can therefore be computed once and reused at every iteration of the V-cycle, mirroring how a classical multigrid method assembles its operator hierarchy once before iterating. The learnable setup phase precomputes and caches these features for all levels, so that each subsequent iteration reduces to a sequence of linear convolutions against the cached coefficients, as in Eq. 5 . This keeps the heavy, nonlinear feature extraction out of the iterative loop and substantially lowers the per-iteration cost. For brevity, we omit the parameter subscript \(\theta\) in the learnable operators below.

3.1.0.1 Feature Hierarchy Construction

We generate two distinct sets of cached physical features for each grid level \(l\):

1. PDE Features (\(\boldsymbol{a}_l\)): These feature maps act as a learned, high-dimensional generalization of finite difference stencil coefficients. They encode local dispersion corrections and wave speeds required to approximate the Helmholtz operator \(A_l\).

2. Smoother Features (\(\boldsymbol{s}_l\)): These maps serve as spatially varying relaxation parameters. They allow the solver to dynamically adjust smoothing intensity based on local heterogeneity.

The detailed construction of the multigrid hierarchy is summarized in Algorithm 1. To facilitate the generation of physical features, we employ a set of non-linear convolutional neural networks (CNNs) at each level. Specifically, let \(f_{\text{lift}}\) denote the lifting network that maps the raw physical coefficients into the high dimensional latent feature space. For level \(l\), we define \(f_{\boldsymbol{s}}^{l}\) as the network projecting PDE features \(\boldsymbol{a}_l\) into smoother features \(\boldsymbol{s}_l\), and \(f_{\boldsymbol{a}}^{l}\) as the restriction network mapping features \(a_l\) to the coarser level. All three types of networks (\(f_{\text{lift}},f_{\boldsymbol{s}}^{l},f_{\boldsymbol{a}}^{l}\)) are composed of multi-channel non-linear convolutional layers. Conversely, the inter-grid transfer operators \(\mathcal{R}_{l}^{l+1}\) and \(\mathcal{P}_{l+1}^{l}\) are strictly multi-channel linear convolutional (or transpose convolutional) layers to preserve the linearity of the solver.

3.2 Learnable Solving Phase↩︎

Building upon the standard Geometric Multigrid framework [34], the solving phase of McMg executes a recursive V-cycle comprising three core operations: neural smoothing, neural restriction, and neural prolongation. These operations reuse the cached physical features precomputed during the learnable setup phase, modulating linear convolutions of the residuals and corrections as in Eq. 5 . We denote the complete neural V-cycle as \(\mathcal{MG}(\boldsymbol{k}, \boldsymbol{r})\), where \(\boldsymbol{k}\) represents the discrete wavenumber, and \(\boldsymbol{r}\) is the input residual. Since the physical coefficients remain fixed throughout the iterative process, we simplify the notation to \(\mathcal{MG}(\boldsymbol{r})\) for brevity.

The overarching objective of the McMg is to function as a high-quality preconditioner that approximates the inverse Helmholtz operator: \[\mathcal{MG}(\boldsymbol{r}) \approx A^{-1} \boldsymbol{r}.\] Crucially, since the inverse operator \(A^{-1}\) is linear, McMg is designed to preserve linearity with respect to \(\boldsymbol{r}\) for fixed \(\boldsymbol{k}\). Nonlinearity enters only through the medium-dependent setup features, not through the residual-to-correction map used during the solve.

3.2.0.1 Neural Smoothing

In contrast to fixed relaxation schemes like Jacobi or Gauss-Seidel, McMg employs a learnable neural smoother built from the cached smoother features \(\boldsymbol{s}_l\). At grid level \(l\) and smoothing iteration \(i\), the neural smoother is defined as: \[\mathcal{S}_{l}^{(i)}(\boldsymbol{x}) = \mathcal{W}^{(i)} \ast \big(\boldsymbol{s}_l \odot (\mathcal{K}_{\boldsymbol{x}} \ast \boldsymbol{x})\big),\] where \(\mathcal{W}^{(i)}\) denotes the kernel of a standard multi-channel linear convolution. Numerical experiments demonstrate that employing iteration-dependent smoother parameters significantly accelerates convergence. Consequently, while the smoother features are shared, we assign a unique set of weights \(\mathcal{W}^{(i)}\) for each smoothing step \(i\). Conversely, to ensure physical consistency across iterations, the wave propagation physics are governed by a single Neural PDE Operator \(\mathcal{A}_{l}\) shared across the entire level: \[\mathcal{A}_{l}(\boldsymbol{x}) = \mathcal{W}_l \ast \boldsymbol{x} + \boldsymbol{a}_l \odot (\mathcal{K}_{\boldsymbol{x}} \ast \boldsymbol{x}),\] where the first term employs a linear convolution with kernel \(\mathcal{W}_l\) to approximate the Laplacian operator, and the second term accounts for the heterogeneous physics via the precomputed PDE features \(\boldsymbol{a}_l\). Combining these operators, the neural smoothing phase is updated via: \[\boldsymbol{e}_l^{(i)} = \boldsymbol{e}_l^{(i-1)} + \mathcal{S}_{l}^{(i)} \left( \boldsymbol{r}_l - \mathcal{A}_{l}(\boldsymbol{e}_l^{(i-1)}) \right).\] Crucially, these operations are performed within a high-dimensional feature space. Unlike the scalar fields used in classical methods, the multi-channel features \((\boldsymbol{e}_l,\boldsymbol{r}_l)\) attach several values to each grid point. The smoother therefore acts not only as a damping step, but also as a medium-dependent channel-mixing operator that can locally modify amplitude, phase, and direction of latent wave packets.

3.2.0.2 Neural Restriction

The objective of the restriction step is to transfer residuals from the fine grid to the coarse grid while preserving the components that are difficult to eliminate locally. While classical geometric multigrid relies on fixed scalar kernels, our neural restriction generalizes this operation to a learnable linear operator implemented as a multi-channel strided convolution, denoted by \(\mathcal{R}_{l}^{l+1}\). This allows the network to project fine-grid latent residual features onto a coarse latent space that can still distinguish phase and directional components. The coarse-level latent residual features \(\boldsymbol{r}_{l+1}\) are computed as: \[\boldsymbol{r}_{l+1} = \mathcal{R}_{l}^{l+1} \left( \boldsymbol{r}_l - \mathcal{A}_{l}(\boldsymbol{e}_l) \right),\]

3.2.0.3 Neural Prolongation

Analogous to the neural restriction, the neural prolongation step employs a learnable linear operator implemented as a multi-channel transposed convolution, denoted by \(\mathcal{P}_{l+1}^{l}\). This operator maps coarse-grid latent error-correction features back to the fine-grid feature space. Unlike standard fixed interpolation, the learned prolongation defines a multi-channel packet dictionary: one coarse spatial node corresponds to several fine-grid oscillatory responses. The fine-grid latent correction features are updated via: \[\boldsymbol{e}_{l} \leftarrow \boldsymbol{e}_{l} + \mathcal{P}_{l+1}^{l} \left( \boldsymbol{e}_{l+1}\right).\]

Integrating these learnable components, we construct the complete McMg V-cycle. This framework propagates multi-channel phase-space information across the grid hierarchy rather than relying on scalar coarse-grid variables. The overall solving procedure is summarized in Algorithm 2. Note that the algorithm employs two linear convolutional modules to interface with the physical field: an encoder \(g_{\text{lift}}\) that maps the initial scalar residual into the high-dimensional latent space, and a decoder \(g_{\text{proj}}\) that reconstructs the final scalar error correction from the latent features. These lifting and projection maps contain no nonlinear activations, so they preserve linearity with respect to the residual input. The overall procedure is illustrated in Fig. 3.

3.3 Training and Inference Strategies↩︎

McMg is trained by minimizing the residual of the discrete Helmholtz system, which we describe below. The same learned operator can also be embedded within the Born series, yielding the Born+McMg variant which follows the neural preconditioning approach of NPBS [36], we defer its derivation to Appendix 6.2.

We consider the discrete Helmholtz system \[A\boldsymbol{u} = \boldsymbol{f},\] where \(A\) denotes the discretized Helmholtz operator and \(\boldsymbol{f}\) is the right-hand side. Given an approximate solution \(\boldsymbol{u}^{(k)}\) with residual \(\boldsymbol{r}_{\mathrm{H}}^{(k)}\), the most direct use of McMg is as a stationary correction scheme that updates the solution directly by the neural correction at each iteration: \[\label{eq:mcmg95iteration} \boldsymbol{r}_{\mathrm{H}}^{(k)} = \boldsymbol{f} - A\boldsymbol{u}^{(k)}, \qquad \boldsymbol{u}^{(k+1)} = \boldsymbol{u}^{(k)} + \mathcal{MG}\big(\boldsymbol{r}_{\mathrm{H}}^{(k)}\big).\tag{6}\] In this setting, McMg itself defines the iterative solver. To train McMg, following [26], [36], we sample random residuals \[\boldsymbol{r} \sim \mathcal{N}(\boldsymbol{0},\boldsymbol{I}),\] where \(\mathcal{N}\) denotes the multivariate standard normal distribution, and minimize the relative residual of the predicted correction: \[\label{eq:loss95fun} \mathcal{L}_{\mathrm{H}} = \frac{1}{B} \sum_{i=1}^{B} \frac{ \left\| A_{(i)} \mathcal{MG}(\boldsymbol{k}_{(i)}, \boldsymbol{r}) - \boldsymbol{r} \right\|_2 }{ \left\| \boldsymbol{r} \right\|_2 },\tag{7}\] where \(B\) is the batch size and \(i\) denotes the sample index. Here, \(A_{(i)}\) is the discrete operator for the \(i\)-th sample, parameterized by the wavenumber \(\boldsymbol{k}_{(i)}\). We retain the explicit dependence \(\mathcal{MG}(\boldsymbol{k}_{(i)}, \boldsymbol{r})\) rather than the simplified form \(\mathcal{MG}(\boldsymbol{r})\), since each sample is associated with a distinct wavenumber. The residuals \(\boldsymbol{r}\) are resampled at each training iteration to expose the model to a broad distribution of residual modes. This objective trains McMg to approximate the inverse of the Helmholtz operator. The same learned operator can also serve as a preconditioner within GMRES.

3.4 The Multi-channel Coarse Space↩︎

The defining feature of McMg, relative to classical scalar multigrid, is that residuals, corrections, and the coarse spaces themselves all reside in a multi-channel latent space rather than on scalar grids. Although wide feature maps are routine in neural networks, here the extra channels serve a specific numerical purpose: they allow a coarse grid to retain local phase-space information that scalar coarsening is forced to discard.

The local phase-space motivation is as follows. In the high-frequency regime, Helmholtz fields and errors are often locally approximated by WKB/geometrical-optics ansatz functions [10], [11] \[e(x)\approx a(x)e^{i\phi(x)},\qquad |\nabla\phi(x)|\approx k(x).\] Thus the unresolved state includes not only the position \(x\) but also the local wave vector \(\xi=\nabla\phi(x)\). A scalar coarse node cannot distinguish wave packets traveling in different directions or with different phases after the spatial grid has been coarsened; this phase-space viewpoint is also consistent with Gaussian-beam and wave-packet descriptions of high-frequency wave propagation [11], [12]. In McMg, the channel index provides a finite-dimensional local representation of these unresolved variables.

Let \(\alpha=(\mathbf{i},c)\) denote a degree of freedom at the coarsest level \(L\), where \(\mathbf{i}\) is a coarse-grid location and \(c\) is a channel index. If \(\boldsymbol{\delta}_{\alpha}\) is the corresponding one-hot tensor, the raw transfer-defined fine-grid response is obtained by the composite learned prolongation \[\label{eq:coarse95basis} \Psi_{\alpha}^{\mathrm{raw}} = g_{\text{proj}}\left(\mathcal{P}_{2}^{1} \circ \mathcal{P}_{3}^{2} \circ \cdots \circ \mathcal{P}_{L}^{L-1} (\boldsymbol{\delta}_{\alpha})\right),\tag{8}\] where \(g_{\text{proj}}\) maps the latent correction back to the physical scalar field. This raw basis furnishes a learned tentative packet dictionary shared across media. The basis effectively used by a V-cycle is further modified by the medium-dependent smoother and neural PDE operators. If \(\mathcal{H}_l(\boldsymbol{k})\) denotes the homogeneous linear part of the level-\(l\) update induced by the cached features \(\boldsymbol{a}_l,\boldsymbol{s}_l\), then a schematic effective response is \[\label{eq:effective95packet95basis} \Psi_{\alpha}^{\mathrm{eff}}(\boldsymbol{k}) \approx g_{\text{proj}}\, \mathcal{H}_1(\boldsymbol{k})\mathcal{P}_2^1 \mathcal{H}_2(\boldsymbol{k})\mathcal{P}_3^2 \cdots \mathcal{P}_L^{L-1}\boldsymbol{\delta}_{\alpha}.\tag{9}\] Thus the transfer operators provide a learned multi-channel packet dictionary, while the cached medium-dependent smoothers and neural PDE operators transform it into an operator-adapted effective coarse representation. This is analogous in spirit to smoothed aggregation AMG, except that the tentative basis is multi-channel and oscillatory rather than scalar and smooth.

For classical piecewise-polynomial bases (e.g., \(P_1\) elements in standard GMG), the coarse mesh size \(H\) is constrained by scalar sampling and dispersion requirements; suppressing the accumulated phase error (the pollution effect) can demand conditions such as \(H k^2 \lesssim C_0\) [1]. Once \(H\) is too large relative to the wavelength, a basis with a single scalar degree of freedom per node can no longer approximate oscillatory error components, and coarse-grid correction becomes ineffective. The standard remedy is to enrich the coarse space: advanced multiscale methods such as Localized Orthogonal Decomposition (LOD) [37]–[41] and the Generalized Multiscale Finite Element Method (GMsFEM) [42]–[44] attach several problem-dependent oscillatory basis functions to each coarse node, which improves the resolution condition.

McMg follows the same enrichment principle, with the per-node basis learned end-to-end rather than constructed analytically. We visualize the learned bases for a six-level model trained on the OpenBreastUS dataset under the configuration of Section 4.1, with one coarse level added beyond the five-level solver used there so as to expose the behavior of the coarsest scale; here \(k_0 D \approx 500\) and the coarsest level \(L=6\) has spacing \(H \approx 32h \approx 5\lambda\). A scalar piecewise-polynomial coarse space would be far too sparse in this regime. As shown in Figure 4, activating a single coarse degree of freedom at level \(L=6\)—a unit impulse in one channel \(c\) at one coarse location \(\mathbf{i}\)—induces a fine-grid response that is not a localized interpolation kernel but a highly oscillatory wave packet spanning roughly \(12\) wavelengths. The coarse representation stays expressive not by reducing \(H\), but by enriching each node with \(C\) learned packet coefficients.

This raises a natural question: are the channels genuinely doing independent work, or is the wide representation largely redundant? To answer it, we fix the center point of the coarsest grid, set one channel at this point to one and all other coarse-level entries to zero, and then map this tensor to the physical grid using Eq. 8 . Repeating this procedure over channels gives a family of fine-grid responses \(\{\Psi^{(c)}\}_{c=0}^{C-1}\). Figure 5 (a) shows representative responses for channels \(0\), \(10\), \(20\), and \(25\), cropped around the heterogeneous scattering region. Although these responses are generated from the same spatial location and the same learned prolongation hierarchy, they exhibit different phase, orientation, and interference patterns. The channels therefore encode distinct physical components rather than replicated scalar bases. We further quantify this diversity by stacking the vectorized responses \(\Psi^{(c)}\) across all channels and computing their singular values. The resulting spectrum, shown in Figure 5 (b), decays gradually rather than collapsing to a few dominant modes. This indicates that the learned channels span a genuinely multi-dimensional coarse subspace. The multi-channel construction is therefore an essential part of the learned coarse representation, rather than a simple increase in feature width.

3.5 Layer-by-Layer Progressive Finetuning for Large Domain Scalability↩︎

Achieving efficient convergence on large-scale computational domains presents a significant challenge for learning-based solvers. Models trained on small domains can accurately approximate local wave interactions but often fail to capture the long-range correlations essential for global error propagation. Consequently, when deployed on larger grids, these locally trained models exhibit degraded convergence rates, necessitating numerous outer iterations. A spectral interpretation of this long-range degradation is provided in Appendix 6.4.

While finetuning the entire model on the target large domain can mitigate this issue, it incurs prohibitive computational costs due to the increased grid resolution and the expense of backpropagating through the full depth of the multigrid hierarchy. To address this trade-off between scalability and efficiency, we propose the Layer-by-Layer Progressive Finetuning (LLPF) strategy. LLPF leverages the hierarchical nature of McMg to systematically extend the spatial support of the learned solver with minimal training overhead.

The LLPF procedure, illustrated in Figure 6, adapts a model trained on a small initial domain \(\mathcal{D}_1\) to a sequence of increasingly larger domains \(\{\mathcal{D}_k\}_{k=2}^K\), where \(\mathcal{D}_K\) represents the target large-scale problem. The process consists of two main phases:

1. Base Model Initialization (\(\mathcal{D}_1\)): We commence by training a standard McMg solver with \(L\) levels on the smallest domain \(\mathcal{D}_1\). This step efficiently establishes the baseline parameters for resolving high-frequency local errors and short-range discrete physics.

2. Recursive Expansion and Finetuning (\(\mathcal{D}_k\)): For each subsequent domain \(\mathcal{D}_k\) larger than \(\mathcal{D}_{k-1}\), we perform the following operations:

   • Architecture Growth: We append an additional, coarser level to the bottom of the V-cycle hierarchy. This extends the effective receptive field of the solver, enabling it to represent the lower-frequency modes characteristic of the larger domain.

   • Progressive Locking: The parameters of all previously trained levels (fine to intermediate grids) are frozen. These levels already capture the local physics correctly; preserving them prevents "catastrophic forgetting" and ensures stability.

   • Selective Optimization: We finetune only the parameters associated with the newly added coarse level using data generated on \(\mathcal{D}_k\). Since the fine-grid errors are already handled by the frozen layers, the optimization focuses exclusively on resolving the residual global error components.

This strategy keeps the training cost tractable as the problem size grows. By restricting backpropagation to a single variable level at each stage, LLPF extends the effective support of the learned Green’s operator without requiring full-model retraining.

4 Numerical Experiments↩︎

In this section, we present a systematic evaluation of the proposed McMg method, benchmarking its performance against established baseline algorithms across diverse heterogeneous datasets. Unless otherwise stated, all computational frameworks were implemented using PyTorch 2.6.0 with CUDA 12.6 acceleration. The experiments were conducted on a high-performance workstation equipped with an Intel Xeon Gold 6444Y CPU and an NVIDIA RTX A6000 GPU.

4.1 Benchmarking McMg against Classical and Neural Solvers↩︎

We evaluate McMg across three challenging regimes, each stressing a different failure mode of classical solvers: a high-frequency setting on the OpenBreastUS dataset [45], a high-contrast setting on the CurveFault-B dataset, and a large-scale three-dimensional setting on the Kimberlina 1.2 CCUS dataset, the latter two from OpenFWI [46]. We deploy the learned McMg both as a standalone solver and as a preconditioner, and compare it against classical and neural solvers. We describe each setting and its training protocol below, then the evaluation protocol, and finally report training cost and solver performance.

4.1.0.1 High-frequency regime (OpenBreastUS).

The OpenBreastUS dataset comprises realistic 2D USCT breast models (Figure 7 (a)) on a \(480\times 480\) grid. We target the challenging low points-per-wavelength (\(\texttt{ppw}\)) regime, where the pollution effect degrades the discrete operator’s spectral properties and destabilizes iterative solvers; here McMg’s multi-channel adaptive operators conform to local heterogeneity and remain stable where fixed-stencil approaches fail. Normalizing the domain to a unit square (\(1\,\text{m}\times 1\,\text{m}\)) gives a grid spacing \(h = 1/480\,\text{m}\) and \(\texttt{ppw} = 2\pi/(k_0 h) = 6\), where \(k_0\) is the background wavenumber. The resulting dimensionless wavenumber \(k_0 D \approx 500\) (\(D\) the domain size) places the problem firmly in the high-frequency regime [6]. We use an absorbing layer of width 48 grid points and multigrid levels \([1,2,4,8,8]\), and train on 1000 samples for 400 epochs with batch size 10 using Adam.

4.1.0.2 High-contrast regime (CurveFault-B).

The CurveFault-B dataset features sharp discontinuities and large velocity variations, from \(c_{\min} = 1500\) m/s to \(c_{\max} = 4500\) m/s (Figure 7 (b)). CBS is challenged here: higher contrasts enlarge the scattering potential and make the constant-background Green’s preconditioner less contractive, leading to slower convergence. McMg instead generates adaptive stencils and smoothers with nonlinear dependence on the local wavenumber \(k(x)\), capturing the wave physics at sharp interfaces and improving conditioning. We normalize the domain to the unit square \([0,1]\times[0,1]\) m on a \(256\times 256\) grid, apply an absorbing layer of width 32 grid points, and set \(\texttt{ppw} = 10\). We train on 3000 samples with batch size 30; all other settings follow the high-frequency configuration.

4.1.0.3 Large-scale 3D regime (Kimberlina 1.2 CCUS).

Three-dimensional problems have long challenged AI-based PDE solvers, whose parameter counts and memory footprints grow rapidly with domain size. Thanks to its Setup and Parameter Sharing designs, McMg both trains and runs inference on \(224\times 224\times 224\) volumes on a single NVIDIA A6000 GPU. We adopt the Kimberlina 1.2 CCUS dataset from OpenFWI, which models subsurface evolution across CO₂ injection years on an original \(400\times 400\times 350\) grid. We interpolate each model to \(192\times 192\times 192\) and add a 16-point absorbing layer, giving a \(224\times 224\times 224\) domain with \(\texttt{ppw}=6\), and sound speeds from \(c_{\min}=2637\) to \(c_{\max}=4201\) m/s—over 11 million degrees of freedom. We train on the first 35 years (630 samples) and validate on the 100th year (50 samples), using batch size 3 and multigrid levels \([1,2,4,6,8]\) with 20 channels; all other settings follow the high-frequency configuration. The resulting model has only 3.82 M parameters and uses roughly 45 GB of GPU memory for training and 6 GB for inference.

4.1.0.4 Baselines and evaluation protocol.

On the high-frequency benchmark, we group our baselines into two classes. The classical solvers are CBS, a GMRES solver preconditioned by shifted-Laplacian multigrid (GMRES+SL) [4], and a GPU-accelerated sparse LU factorization as a direct-solver reference; the neural solver is MGCFNN [26]. We evaluate McMg as a standalone solver, as a GMRES preconditioner (GMRES+McMg), and as a Born-series preconditioner (Born+McMg). To ensure a fair comparison, all solvers act on the same five-point finite-difference discretization; the Born-series solvers (CBS and Born+McMg) retain their FFT-based Green operator but construct it from the discrete five-point symbol rather than the continuous Laplacian symbol, matching the discretization of the others (Appendix 6.2.2). All iterative methods except CBS run in mixed precision, applying the preconditioner in FP32 and evaluating the residual in FP64; CBS runs entirely in FP64, as its preconditioning reduces to a single element-wise multiplication for which precision conversion would cost more than it saves. All runs share the same GPU environment, with iteration counts and wall-clock times averaged over 50 test samples at a relative tolerance of \(1\times10^{-6}\). MGCFNN uses the same configuration as McMg. The sparse LU reference is implemented with the highly optimized NVIDIA cuDSS¹ library in single precision, and the convolutional forward pass is accelerated with PyTorch’s just-in-time (JIT) compilation. The full parameterization of the classical baselines—the shifted-Laplacian shift, smoother, coarse solver, and GMRES restart length, together with the CBS background wavenumber and relaxation—is recorded in Appendix 6.3.

4.1.0.5 Training cost.

Before turning to solver performance, we compare the training cost of McMg and MGCFNN on the high-frequency benchmark under an identical protocol (Table 1). McMg attains lower training and validation losses while using roughly an eighth of the parameters and a third of the per-epoch time, reflecting the efficiency of its parameter-sharing multigrid design.

4.1.0.6 Convergence and runtime.

Table ¿tbl:tab:cbs95mpmg95precision? reports solver performance across the three regimes; we begin with the high-frequency benchmark. Among the classical methods, CBS attains the lowest wall-clock time (\(0.372\) s), ahead of shifted-Laplacian-preconditioned GMRES (GMRES+SL, \(3.290\) s) and a highly optimized sparse LU direct solver (\(1.747\) s); we therefore adopt it as the reference classical solver throughout. Relative to CBS, McMg reduces the iteration count by \(110.6\times\); because a single McMg iteration is more expensive than a single CBS iteration, this translates into a \(10.3\times\) reduction in wall-clock time. That this gain originates in the learned preconditioner is confirmed by a controlled comparison: Born+McMg uses the same Born iteration as CBS and differs only in its preconditioner, yet converges in \(8.0\) iterations as compared with \(929.2\) for CBS. McMg likewise compares favorably with the neural baseline, converging nearly \(4\times\) faster than MGCFNN. Interestingly, using McMg as a GMRES preconditioner reduces the iteration count further, from \(8.4\) to \(7.2\), but does not improve the wall-clock time (\(0.036\) s versus \(0.059\) s), indicating that the standalone McMg architecture is already efficient enough that the overhead of the Krylov recurrence outweighs the benefit of fewer iterations. The same advantage extends to the remaining regimes: relative to CBS, McMg reduces the iteration count by \(51.3\times\) and \(45.1\times\) and the wall-clock time by \(3.4\times\) and \(6.8\times\) in the high-contrast and three-dimensional settings, respectively, confirming its effectiveness across high-frequency, high-contrast, and large-scale regimes. Figures 8, 9, and 10 corroborate these findings, exhibiting accurate wavefields and rapid, monotone residual decay.

4.2 Domain Generalization via Patch-Based Training↩︎

The generalization capability of McMg is underpinned by the intrinsic locality of the learned stencils. Since the stencil coefficients at a given point are primarily governed by the local wave velocity structure and discretization errors rather than by the total number of grid points, the learned map can transfer from small patches to larger domains. This property enables a train-on-patches, deploy-on-full-domains strategy, where McMg is trained on computationally inexpensive subdomains and then tested on larger domains without finetuning.

To validate this capability, we train McMg exclusively on small sub-patches (\(256 \times 256\)) randomly extracted from the Marmousi velocity model (Figure 11a). Subsequently, the trained model is employed to solve the Helmholtz equation on the complete \(1024 \times 1024\) domain (Figure 11b). The experimental setup follows the configuration in Section 4.1.0.2, utilizing a sampling density of \(\texttt{ppw} = 12\). The equation is discretized using a standard 5-point FDM, augmented with a sponge layer of width 32 to suppress boundary reflections.

Figure 12 and Figure 13 illustrate the solver’s convergence performance on a cropped validation patch and the full domain, respectively. As anticipated, the model demonstrates robust generalization: the number of iterations required for the full-scale problem exhibits only a marginal increase compared to the sub-problems. This confirms that the local operators learned by McMg effectively capture the underlying physics, regardless of the global scale of the simulation.

4.3 Scalability and Cross-Scale Generalization↩︎

A pivotal advantage of our McMg is the robust generalization across scales and wavenumbers. Our model, trained on small-scale domains, can be directly applied to large-scale problems with significantly higher wavenumbers without retraining. This “train on small-scale, deploy on large-scale” capability represents a substantial advancement over neural operators that typically require retraining when the problem size or frequency range changes. We investigate this scalability by training McMg on a limited-scale dataset (\(480\times480\)) and subsequently deploying it on problems of exponentially increasing size and frequency, ranging from a grid size of \(512 \times 512\) (with \(k_0 D \approx 500\)) up to \(8192 \times 8192\) (with \(k_0 D \approx 8000\)). Based on these results, we further improve large-domain scalability using the LLPF strategy described in Section 3.5. The experiments are conducted on the OpenBreastUS Breast dataset, utilizing the same configuration detailed in Section 4.1 with a convergence tolerance of \(10^{-6}\).

To evaluate scalability, we compare McMg, Born+McMg, and the CBS solver (Table 2). For McMg, the per-iteration cost scales linearly as \(\mathcal{O}(N)\) with respect to the total number of grid points \(N\). For Born+McMg and CBS, which apply the Helmholtz operator via the FFT, the per-iteration cost is \(\mathcal{O}(N\log N)\): the \(\log N\) factor stems solely from the FFT-based operator application, while all remaining operations scale linearly in \(N\).

The iteration count nonetheless increases monotonically with the domain size. To see why, recall the McMg correction step Eq. 6 , so the correction applied at step \(k\) is \(\Delta\boldsymbol{u}^{(k)} = \mathcal{MG}(\boldsymbol{r}_{\mathrm{H}}^{(k)})\). Driving the solver with a point source \(\boldsymbol{f} = \boldsymbol{\delta}_{x_s}\) from a zero initial guess, the first correction \[\label{eq:green95visual} \Delta\boldsymbol{u}^{(0)} = \mathcal{MG}(\boldsymbol{\delta}_{x_s}) \approx A^{-1}\boldsymbol{\delta}_{x_s}\tag{10}\] is exactly the learned approximation of the discrete Green’s function centered at the source \(x_s\), and each subsequent correction applies the same operator to the remaining residual, transporting the wave further across the domain. Visualizing the first five corrections \(\{\Delta\boldsymbol{u}^{(k)}\}_{k=0}^{4}\) (Figure 14) therefore exposes the effective spatial support of the learned Green’s function. A model trained only on a small domain reproduces a spatially truncated Green’s function (Figure 14 (a)): it acts as a localized operator that captures only short-range interactions, so each correction reaches a limited range and many outer iterations are needed to propagate the solution across a larger domain. To mitigate this growth in iteration count, we next consider finetuning and the LLPF strategy.

We evaluate different finetuning strategies, with results reported in Table ¿tbl:tab:mpmg95bs95scalability95comp95finetune95and95llpf?. Direct finetuning on the \(1024 \times 1024\) domain with \(7\) levels yields the best overall performance, balancing optimization stability and effective Green’s function support. In contrast, finetuning directly on the \(2048 \times 2048\) domain often degrades convergence due to increased optimization difficulty (see Appendix 6.4). LLPF mitigates this issue by freezing existing levels and finetuning only the newly added coarsest level. As shown in Table ¿tbl:tab:mpmg95bs95scalability95comp95finetune95and95llpf?, LLPF consistently improves convergence as coarse levels are added, matching the convergence of direct finetuning on the \(1024 \times 1024\) domain at a small fraction of its cost. Because backpropagation is confined to the single newly added level, extending a five-level model to six levels with LLPF takes only \(0.88\)h, about \(3\times\) faster than directly finetuning a full six-level model on the same \(1024\times1024\) domain (\(2.51\)h), with a comparable reduction in memory. The effective support of the learned Green’s function tells the same story (Figure 14): although LLPF updates only the added level, its correction (Figure 14 (c)) attains the same long-range support as direct finetuning of the full model on the \(1024\times1024\) domain (Figure 14 (b)), and both extend markedly beyond the short-range, truncated correction of the model pretrained on \(512\times512\) (Figure 14 (a)).

4.4 Spectral Convergence Analysis↩︎

We analyze the error evolution in the frequency domain to investigate how McMg handles spectral error modes. A major limitation of GMG is its inability to attenuate characteristic components [25] (error modes with wavenumber \(|\mathbf{\xi}| \approx k\)). These modes are difficult to eliminate because they generate small residuals for local smoothers and suffer from severe phase errors on scalar coarse grids.

We examine a test sample with a heterogeneous breast sound speed map (Figure 15 (a)). Figures 15 (b) and 15 (e) visualize the spectral error magnitude at iterations 1 and 7, respectively. The initial error (Figure 15 (b)) is heavily concentrated around the characteristic ring \(|\mathbf{\xi}| \approx k\), confirming that these modes dominate the error landscape. McMg suppresses these components by iteration 7 (Figure 15 (e)), while also reducing the lower- and higher-frequency bands.

The quantitative breakdown in Figure 15 (f) confirms this observation. The error at the characteristic frequency \(k\) decays at an approximately exponential rate, parallel to both the low-frequency (\(k/10\)) and high-frequency (\(2k\)) components. This comparable decay across diverse bands supports the interpretation that McMg’s multi-channel coarse representation and adaptive operators preserve the wave information needed to reduce characteristic modes that typically stagnate standard GMG solvers.

5 Conclusions↩︎

We have presented McMg, a learned multigrid preconditioner for heterogeneous Helmholtz equations that replaces end-to-end solution prediction with iterative residual-to-correction learning. The main architectural idea is to replace scalar coarse-grid correction by a multi-channel phase-space coarse representation: physical space is coarsened, but each coarse node carries latent wave-packet coefficients that can represent local phase, direction, and scattering information. Medium-dependent neural PDE operators and smoothers then adapt this packet representation to the local wave speed. Two design choices make the method practical at scale. First, parameter sharing within the multigrid hierarchy yields a compact model. Second, the setup phase precomputes and caches the nonlinear physical features once per medium, so that each iteration reduces to a sequence of linear convolutions with cached coefficients. McMg can also be incorporated into the CBS framework, yielding Born+McMg. This variant is particularly effective in strongly scattering media, where training in the preconditioned Born-series residual space improves stability and convergence efficiency. To extend the solver to large domains without full retraining, LLPF enlarges the effective support of the learned Green’s operator by adding and finetuning only new coarse levels.

Across high-frequency (OpenBreastUS), high-contrast (OpenFWI CurveFault-B), and large-scale three-dimensional (Kimberlina) benchmarks, McMg converges in substantially fewer iterations and shorter wall-clock time than the fastest traditional solver in our tests, and outperforms recent neural solvers at comparable accuracy. Relative to CBS, this strongest traditional reference, McMg reduces iteration counts by \(110.6\times\), \(51.3\times\), and \(45.1\times\) on the three benchmarks, with corresponding wall-clock speedups of \(10.3\times\), \(3.4\times\), and \(6.8\times\). Spectral and coarse-space analyses further indicate that the learned operators reduce characteristic error modes near \(|\boldsymbol{\xi}|\approx k\) and, by attaching a multi-channel local basis to each coarse node, represent oscillatory components on grids where a scalar coarse space would be too sparse.

Several directions remain open. A local Fourier or Wigner diagnostic should directly test whether the learned channels tile the characteristic shell \(|\xi|\approx k(x)\); effective coarse-space approximation tests should compare McMg against scalar GMG, LOD, and GMsFEM bases; and stability diagnostics should report field-of-values, pseudospectral, or worst-case residual amplification measures for the learned preconditioner. The present formulation relies on structured grids, so extending the learned stencils to unstructured meshes and complex geometries is important. Finally, McMg should be embedded as an inner solver in large-scale inverse problems, where repeated Helmholtz solves for many sources and frequencies are a dominant computational cost.

6 Additional Derivations and Experiments↩︎

6.1 Alternative Discretizations and Boundary Conditions↩︎

The main text discretizes the Helmholtz operator with the five-point finite-difference (FDM) stencil Eq. 3 and truncates the domain with a sponge layer. This appendix details the two principal alternatives, the perfectly matched layer (PML) as an absorbing boundary and the Fourier spectral method (FSM) as a discretization, and then compares all four combinations empirically.

6.1.1 Derivation of the Perfectly Matched Layer (PML)↩︎

The Perfectly Matched Layer (PML) method extends the spatial coordinates into the complex plane to introduce artificial decay without reflection. For the Helmholtz equation, this is achieved by applying a complex coordinate stretching: \[\tilde{x}_j = x_j + i p_j(x_j), \quad j = 1, \dots, D,\] where \(p_j(x_j)\) is a monotonically increasing function in the \(j\)-th coordinate direction. The corresponding Jacobian of the transformation is defined as \[\frac{\partial \tilde{x}_j}{\partial x_j} = 1 + i \frac{\partial p_j(x_j)}{\partial x_j}, \quad j = 1, \dots, D.\] Introducing the complex stretching factor: \[s_j(x_j) = 1 + i \frac{\partial p_j(x_j)}{\partial x_j} = 1 + i\sigma_j(x_j),\] where \(\sigma_j(x_j)\) denotes the damping profile in the \(j\)-th direction, we have \(\sigma_j(x_j) = 0\) inside the physical domain \(\Omega\) and \(\sigma_j(x_j) > 0\) within the absorbing layer \(\Omega_{\text{abl}}\) (often chosen as a smooth polynomial). The derivative with respect to the complex coordinate becomes \[\frac{\partial}{\partial \tilde{x}_j} = \frac{1}{s_j(x_j)} \frac{\partial}{\partial x_j}, \quad j = 1, \dots, D.\] Consequently, the Laplacian operator modified by the PML becomes \[\tilde{\Delta} = \sum_{j=1}^D \frac{1}{s_j(x_j)} \frac{\partial}{\partial x_j} \left( \frac{1}{s_j(x_j)} \frac{\partial}{\partial x_j} \right) = \sum_{j=1}^D \left( \frac{1}{s_j^2(x_j)} \frac{\partial^2}{\partial x_j^2} - \frac{s_j'(x_j)}{s_j^3(x_j)} \frac{\partial}{\partial x_j} \right).\] The Helmholtz equation with PML is then written as \[-\tilde{\Delta} u(x) - k(x)^2 u(x) = f(x), \quad x \in \Omega \cup \Omega_{\text{abl}}.\]

6.1.2 Fourier Spectral Discretization↩︎

As an alternative to the finite difference discretization used in the main text (Eq. 2 ), the Fourier spectral method (FSM) discretizes the Helmholtz operator globally by leveraging the spectral representation of the Laplacian: \[\Delta u(x) = -\mathcal{F}^{-1} \!\left[\, |\xi|^2 \, \mathcal{F}[u(x)] \,\right],\] where \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) denote the Fourier and inverse Fourier transforms, respectively, and \({\xi}\in \mathbb{R}^D\) is the Fourier wavenumber variable. Applying this to the Helmholtz equation yields a formulation that is exact for the Laplacian term under periodic boundary conditions: \[\mathcal{F}^{-1} \left[ |\xi|^2 \, \mathcal{F}[u(x)] \right] - k(x)^2 u(x) = f(x), \quad x \in \mathbb{R}^{D}.\] Discretization via the Fast Fourier Transform (FFT) yields: \[\label{eq:FSM95discrete} \mathrm{FFT}^{-1} \left[ |\boldsymbol{\xi}|^2 \, \mathrm{FFT}[\boldsymbol{u}] \right] - \boldsymbol{k}^2 \boldsymbol{u} = \boldsymbol{f}.\tag{11}\] Unlike the finite difference method, FSM minimizes dispersion errors by resolving derivatives globally. However, the FFT inherently imposes periodic boundary conditions, which can introduce wrap-around artifacts in non-periodic scattering problems unless combined with sufficient damping (PML).

6.1.3 Empirical comparison of discretizations and boundaries↩︎

We compare the four combinations of discretization (FDM, FSM) and absorbing boundary (sponge, PML) under the model and dataset configuration of Section 4.1.0.1, with results in Table 3. Two trends emerge. Within each discretization, the PML degrades both training loss and convergence relative to the sponge layer: it modifies the differential operator itself, replacing the Laplacian by the complex-stretched operator \(\tilde{\Delta}\) derived above, whose additional first-order terms make the learned operator harder to optimize, whereas the sponge layer acts only on the zeroth-order term and leaves the Laplacian unchanged. Across discretizations, FSM is consistently less favorable than FDM under the same boundary, since the periodic boundaries implied by the FFT admit wrap-around error when the attenuation is not sufficiently strong. The two effects compound, making FSM+PML by far the hardest configuration, while FDM with a sponge layer, the setting used throughout the main text, is the most efficient.

6.2 The Born-Series Variant and Its Discrete Green Operator↩︎

This appendix collects two technical components referenced from the main text: the Born-series formulation of McMg (the Born+McMg variant of Section [sec:mpmg95iteration95scheme]) and the discrete Green operator used to keep all solvers on a common finite-difference discretization in Section 4.1.

6.2.1 McMg for the Born-series system↩︎

For ill-conditioned Helmholtz problems, training McMg directly in the original residual space can lead to unstable optimization and slow convergence. Neural Preconditioned Born Series (NPBS) [36] introduced a unified framework for training neural preconditioners in the preconditioned Born-series residual space, improving both training stability and solver convergence for ill-conditioned Helmholtz systems.

To define this formulation, let \(\boldsymbol{G}\) denote the convolution operator with the background Green’s function associated with the constant background wavenumber \(k_0\), i.e., \(\boldsymbol{G} = (-\Delta - k_0^2)^{-1}\), and let \(\boldsymbol{V}\) be the scattering potential, \[\boldsymbol{V}(\boldsymbol{x}) = k(\boldsymbol{x})^2 - k_0^2.\] The corresponding Lippmann–Schwinger equation [47], [48] is \[\boldsymbol{u} = \boldsymbol{G}(\boldsymbol{V}\boldsymbol{u}+\boldsymbol{f}).\] Equivalently, \[(\boldsymbol{I}-\boldsymbol{G}\boldsymbol{V})\boldsymbol{u} = \boldsymbol{G}\boldsymbol{f}, \qquad \boldsymbol{I}-\boldsymbol{G}\boldsymbol{V}=\boldsymbol{G} A.\] Iterating this fixed-point equation gives the classical Born iteration. Given an approximate solution \(\boldsymbol{u}^{(k)}\), the Born-series residual is \[\boldsymbol{r}_{\mathrm{BS}}^{(k)} = \boldsymbol{G}(\boldsymbol{V}\boldsymbol{u}^{(k)}+\boldsymbol{f}) - \boldsymbol{u}^{(k)} = \boldsymbol{G}(\boldsymbol{f}-A\boldsymbol{u}^{(k)}).\] Thus the Born-series residual is precisely the Helmholtz residual after left preconditioning by \(\boldsymbol{G}\). The classical Born iteration may diverge under strong scattering. The Convergent Born Series (CBS) [7] restores convergence by introducing a pointwise (diagonal) relaxation operator \(\boldsymbol{\gamma}\): \[\label{eq:cbs95eq} \boldsymbol{u}^{(k+1)} = \boldsymbol{u}^{(k)} + \boldsymbol{\gamma} \left( \boldsymbol{G}(\boldsymbol{V}\boldsymbol{u}^{(k)}+\boldsymbol{f}) - \boldsymbol{u}^{(k)} \right).\tag{12}\] With a suitable choice of \(\boldsymbol{\gamma}\), this iteration is guaranteed to converge under the CBS framework, although convergence can still be slow in high-contrast media. Following NPBS, we replace this fixed diagonal relaxation with a learned McMg: \[\boldsymbol{u}^{(k+1)} = \boldsymbol{u}^{(k)} + \mathcal{MG} \left( \boldsymbol{G}(\boldsymbol{V}\boldsymbol{u}^{(k)}+\boldsymbol{f}) - \boldsymbol{u}^{(k)} \right).\] The training objective is chosen to match the preconditioned residual space. We sample residuals \(\boldsymbol{r}\) as in Section [sec:mpmg95iteration95scheme] and use \(\boldsymbol{G}_{(i)}\boldsymbol{r}\) as the network input, so that training and inference use the same preconditioned coordinates. The corresponding loss is \[\mathcal{L}_{\mathrm{BS}} = \frac{1}{B} \sum_{i=1}^{B} \frac{ \left\| \left(\boldsymbol{I}-\boldsymbol{G}_{(i)}\boldsymbol{V}_{(i)}\right) \mathcal{MG} \left( \boldsymbol{k}_{(i)}, \boldsymbol{G}_{(i)}\boldsymbol{r} \right) - \boldsymbol{G}_{(i)}\boldsymbol{r} \right\|_2 }{ \left\| \boldsymbol{G}_{(i)}\boldsymbol{r} \right\|_2 }.\] This trains McMg to approximate the inverse of the preconditioned operator \(\boldsymbol{G} A\), rather than the inverse of the original Helmholtz operator \(A\). The left preconditioning by \(\boldsymbol{G}\) regularizes the spectrum of the operator to be inverted, which is what stabilizes training and accelerates convergence in the ill-conditioned regime.

6.2.2 Discrete Green operator on the finite-difference grid↩︎

The Born-series solvers apply \(\boldsymbol{G}\) through the FFT, which makes them especially efficient but, in their original form, ties them to the spectral symbol \(-|\boldsymbol{\xi}|^2\) of the continuous Laplacian. To compare CBS and Born+McMg against the finite-difference solvers on an identical discrete system (Section 4.1), we instead build \(\boldsymbol{G}\) from the symbol of the discrete Laplacian, retaining the FFT purely as a fast diagonalization tool rather than as a spectral discretization. For the standard five-point stencil in two dimensions with grid spacing \(h\), the discrete Laplacian is diagonalized by the Fourier modes with eigenvalues \[\label{eq:discrete95symbol} \lambda_h(\xi_x,\xi_y) = -\frac{4}{h^2}\sin^2\!\left(\frac{\xi_x h}{2}\right) -\frac{4}{h^2}\sin^2\!\left(\frac{\xi_y h}{2}\right),\tag{13}\] which is the second-order consistent approximation of \(-|\boldsymbol{\xi}|^2\) and recovers it as \(h\to0\). Replacing \(-|\boldsymbol{\xi}|^2\) by \(\lambda_h\) in the background Green operator, \(\boldsymbol{G} = (-\Delta_h - k_0^2)^{-1}\) is realized as multiplication by \((-\lambda_h - k_0^2)^{-1}\) in Fourier space, so the resulting Green’s function is consistent with the five-point finite-difference operator. Consequently, the LU factorization, the shifted-Laplacian preconditioner, and the Born-series methods all act on the same discrete Helmholtz system, and the comparison reflects differences between solvers rather than between discretizations.

6.2.3 Native continuous-spectrum CBS and Born+McMg↩︎

The main-text comparison (Section 4.1) deliberately keeps every solver on the common five-point discretization through the discrete Green operator above. For completeness, we also evaluate the Born-series solvers in their native form, where \(\boldsymbol{G}\) is built from the continuous Laplacian symbol \(-|\boldsymbol{\xi}|^2\) and applied via the FFT, rather than from the discrete symbol \(\lambda_h\). We retrain Born+McMg in this setting and compare it against native CBS under the OpenBreastUS high-frequency configuration of Section 4.1.0.1, so that the only difference from CBS is the preconditioner. As shown in Table 4, native Born+McMg converges in \(6.4\) iterations against \(772.1\) for native CBS, a \(120.6\times\) reduction in iteration count and a \(12.0\times\) wall-clock speedup. These gains match, and slightly exceed, those observed on the common FDM discretization in the main text (Table ¿tbl:tab:cbs95mpmg95precision?), confirming that the advantage of the learned preconditioner persists in CBS’s native continuous-spectrum setting.

6.3 Baseline Solver Configurations↩︎

For reproducibility, we record the parameters of the classical baselines compared in Section 4.1. All baselines act on the same five-point finite-difference discretization (Appendix 6.2.2) and run in the same GPU environment; the precision of each method is chosen as described in Section 4.1, with CBS in FP64 and the remaining iterative methods in mixed precision. Hardware is therefore identical across methods, while the precision of each solver is set to its own fastest configuration so that none is disadvantaged.

6.3.0.1 Shifted-Laplacian-preconditioned GMRES (GMRES+SL).

The preconditioner is the complex-shifted Laplacian operator [4] \[\label{eq:csl95shift} M_{(\beta_1,\beta_2)} = -\Delta_h - (\beta_1 + i\,\beta_2)\,k^2,\tag{14}\] with shift \((\beta_1,\beta_2) = (1,\,1.0)\), i.e.\(M = -\Delta_h - (1+i)\,k^2\). The shifted system is approximately inverted by a single multigrid V-cycle whose smoother is weighted Jacobi with relaxation \(\omega = 0.66\) applied to the diagonal of \(M\), using \((\nu_1,\nu_2) = (2,2)\) pre- and post-smoothing steps; the coarse level is solved by the same weighted-Jacobi iteration for \(2\) steps rather than by a direct solve. The comparatively strong imaginary shift \(\beta_2 = 1.0\) (versus the textbook value \(0.5\)) is adopted precisely because the coarse problem is handled by this inexpensive inner iteration: the added dissipation keeps the V-cycle stable and contractive without an exact coarse solve. GMRES is restarted every \(25\) iterations, i.e.GMRES(\(25\)), which bounds the Krylov-subspace memory. We note that GMRES+McMg converges in \(7.2\) iterations (Table ¿tbl:tab:cbs95mpmg95precision?), well below the restart length, so restarting is never triggered for McMg and influences only the SL baseline.

6.3.0.2 Convergent Born Series (CBS).

CBS is evaluated using the standard aggressive parameterization of Osnabrugge et al. [7]. With the scattering potential \(\boldsymbol{V}(\boldsymbol{x}) = k(\boldsymbol{x})^2 - k_0^2\) of Appendix 6.2, we center the background wavenumber on the real part of the squared-wavenumber range, \(k_0^2 = \tfrac12\big(\min_{\boldsymbol{x}} \mathrm{Re}\{k(\boldsymbol{x})^2\} + \max_{\boldsymbol{x}} \mathrm{Re}\{k(\boldsymbol{x})^2\}\big)\), the standard choice that reduces the required absorption; for the real-valued \(k^2\) of our benchmarks this minimizes \(\max_{\boldsymbol{x}}|\boldsymbol{V}(\boldsymbol{x})|\) exactly. CBS then introduces a uniform absorption \(\epsilon\), using the complex potential \(\boldsymbol{V}_\epsilon(\boldsymbol{x}) = \boldsymbol{V}(\boldsymbol{x}) - i\epsilon\) and the pointwise diagonal preconditioner \[\label{eq:cbs95gamma} \boldsymbol{\gamma}(\boldsymbol{x}) = \frac{i\,\boldsymbol{V}_\epsilon(\boldsymbol{x})}{\epsilon} = 1 + \frac{i\,\boldsymbol{V}(\boldsymbol{x})}{\epsilon};\tag{15}\] the same \(i\epsilon\) is absorbed into the background Green operator, so the Helmholtz operator \(A\) is left unchanged. We take \(\epsilon = \max_{\boldsymbol{x}}|\boldsymbol{V}(\boldsymbol{x})|\), the smallest value permitted by the sufficient CBS convergence condition \(\epsilon \ge \max_{\boldsymbol{x}}|k(\boldsymbol{x})^2 - k_0^2|\) for the no-gain media and absorbing boundaries considered here. The smaller \(\epsilon\) is, the longer the effective range of the absorbing background Green’s function, and the CBS pseudo-propagation distance per iteration scales as \(2k_0/\epsilon\); hence a smaller admissible \(\epsilon\) generally leads to faster convergence. This is therefore the most aggressive choice admitted by the CBS formulation, placing CBS in a favorable standard setting rather than at a parameterization proven to be globally optimal.

6.4 Spectral Analysis of Large Domain Scalability↩︎

Consider the constant-coefficient Helmholtz equation \(-\Delta u - k^2 u = f\), discretized as the linear system \(A \boldsymbol{u} = \boldsymbol{f}\). We aim to learn a parameterized linear operator \(\mathcal{MG} \approx A^{-1}\) by minimizing the unsupervised residual objective \[\min \; \mathcal{L} = \mathbb{E}_{\boldsymbol{r} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{I})} \big\| A \mathcal{MG}(\boldsymbol{r}) - \boldsymbol{r} \big\|_2^2. \label{eq:residual95training95obj}\tag{16}\] To expose the spectral weighting, consider the idealized case in which \(A\) is diagonalized by an orthonormal Fourier or sine basis \(\{\boldsymbol{\phi}_i\}_{i=1}^N\), \[A\boldsymbol{\phi}_i = \lambda_i \boldsymbol{\phi}_i, \qquad \lambda_i \simeq |\boldsymbol{\xi}_i|^2-k^2,\] where the last relation denotes the corresponding discrete Helmholtz symbol. Boundary and PML effects are suppressed in this local Fourier argument. We further consider the diagonal Fourier response of the learned operator, \[\mathcal{MG}(\boldsymbol{\phi}_i) = m_i \boldsymbol{\phi}_i,\] which is exact for a shift-invariant linear operator and serves as the standard local approximation for a trained convolutional multilevel operator. For an input \(\boldsymbol{r}=\sum_i \eta_i \boldsymbol{\phi}_i\), with standardized white-noise coefficients satisfying \(\mathbb{E}[\eta_i\overline{\eta_j}]=\delta_{ij}\), the residual error becomes \[A \mathcal{MG}(\boldsymbol{r}) - \boldsymbol{r} = \sum_{i=1}^N \eta_i \left( \lambda_i m_i - 1 \right)\boldsymbol{\phi}_i .\] Taking the expectation with respect to \(\boldsymbol{r}\) and exploiting the orthogonality of eigenmodes yields: \[\begin{align} \mathcal{L} &= \mathbb{E} \left\| \sum_{i=1}^N \eta_i \left( \lambda_i m_i - 1 \right)\boldsymbol{\phi}_i \right\|_2^2 \\ &= \mathbb{E} \left[ \sum_{i=1}^N \sum_{j=1}^N \eta_i \overline{\eta_j} \left(\lambda_i m_i - 1\right) \overline{\left(\lambda_j m_j - 1\right)} \left\langle \boldsymbol{\phi}_i, \boldsymbol{\phi}_j \right\rangle \right] \\ &= \sum_{i=1}^N \sum_{j=1}^N \mathbb{E}[\eta_i \overline{\eta_j}] \left(\lambda_i m_i - 1\right) \overline{\left(\lambda_j m_j - 1\right)} \left\langle \boldsymbol{\phi}_i, \boldsymbol{\phi}_j \right\rangle \\ &= \sum_{i=1}^N \left| \lambda_i m_i - 1 \right|^2 \\ &= \sum_{i=1}^N |\lambda_i|^2 \left| m_i - \lambda_i^{-1} \right|^2, \end{align} \label{eq:spectral95loss}\tag{17}\] for nonzero \(\lambda_i\). Thus the residual objective controls the residual error \(|\lambda_i m_i-1|\), while the same objective weights the absolute error in the learned inverse response, \(|m_i-\lambda_i^{-1}|\), by \(|\lambda_i|^2\). For the Helmholtz operator, near-resonant or near-propagating modes satisfy \(|\boldsymbol{\xi}_i|^2 \approx k^2\) and hence \(|\lambda_i| \ll 1\) in this idealized representation. These modes have large inverse response and are closely associated with the long-range components of the Green’s function[49], [50]. Consequently, finite-capacity and finite-time training can fit the residual objective while still leaving appreciable errors in the inverse action on these near-resonant components. The point is not that the exact minimizer of the residual loss differs from \(A^{-1}\); rather, the residual formulation gives weak leverage on the large-gain components that are most important for global wave transport.

On a finite domain of length \(L\), the discrete frequency spectrum has resolution \(\Delta \xi = \mathcal{O}(2\pi/L)\). On a small domain, this coarse resolution leaves relatively few near-resonant modes, and the resonance gap \(\delta_N=\min_i |\lambda_i|\) is typically larger. Training on such domains is therefore better conditioned and can learn accurate local correction behavior. However, the learned operator is exposed only to wave interactions over the physical scale of the training domain; when transferred to a larger domain, it behaves like a localized or spatially truncated approximation of the Green’s function and cannot directly correct all long-range error components. By contrast, as \(L\) increases, the spectrum becomes denser near the propagating shell, the number of near-resonant modes grows, and \(\delta_N\) decreases. Direct training on the large domain must then learn many large-gain inverse responses that are weakly weighted in the residual-induced inverse-error metric. In practice, this makes optimization less stable and can lead to a learned operator whose residual loss is small but whose global wave coupling remains inaccurate. The resulting long-range error accumulation explains the degradation in convergence observed when the learned correction is used as an iterative solver, and motivates LLPF as a way to expand the learned Green’s-function support through newly added coarse levels while preserving the local correction already learned on smaller domains.

6.5 Ablation Study on Parameter Sharing↩︎

To assess the impact of our architectural design choices on computational efficiency and model capacity, we conduct an ablation study comparing the proposed McMg against two variants with reduced parameter sharing. The variants are defined as follows:

• McMg: Adopts the full parameter sharing strategy. The Neural PDE Operator \(\mathcal{A}_{l,\boldsymbol{\theta}}\) is shared across all smoothing steps within a level, and physical feature maps (\(\boldsymbol{a}_l, \boldsymbol{s}_l\)) are pre-computed and cached (Setup Phase).

• McMg-1: Disables the sharing of the Neural PDE Operator \(\mathcal{A}_{l,\boldsymbol{\theta}}\). Independent operators are instantiated for every smoothing iteration, though the physical feature maps remain shared.

• McMg-2: Further disables the sharing of physical feature maps (\(\boldsymbol{a}_l, \boldsymbol{s}_l\)). This effectively removes the cached Setup Phase, requiring the network to re-compute physical features from the wavenumber \(\boldsymbol{k}\) at every iteration.

The training characteristics are summarized in Table 5. As expected, reducing parameter sharing (McMg-1 and McMg-2) significantly increases the model size, with McMg-2 requiring nearly \(3\times\) the parameters of the proposed McMg. While McMg-2 achieves a marginally lower training loss due to its increased expressive capacity, this comes at the cost of a substantially longer training time per epoch (35.22s vs. 21.04s).

Table 6 details the inference performance, averaged over 50 test samples with a stopping tolerance of \(1 \times 10^{-6}\). The results highlight a critical trade-off: while the highly parameterized McMg-2 requires slightly fewer iterations to converge (7.5 vs. 8.1), its wall-clock inference time is nearly double that of the proposed method (0.091s vs. 0.049s). This confirms that the computational overhead of re-evaluating physical features and applying unshared operators outweighs the marginal gain in convergence speed. The proposed McMg strikes the optimal balance, delivering comparable convergence with significantly higher computational efficiency.

6.6 Ablation study on the number of multigrid levels↩︎

We investigate the effect of varying the number of multigrid levels in the McMg. To ensure a fair comparison, the total number of smoothing steps is fixed across all levels. The detailed configurations are summarized in Table 7. For instance, the setting \([1,16]\) corresponds to a two-level model, where level \(1\) performs \(1\) smoothing step and level \(2\) performs \(16\) steps. To avoid the influence of other factors, the number of channels at each level is fixed to \(20\). The wavenumber of the Helmholtz equation and the dataset follow the same settings as in Section 4.1.0.1, while here the discretization is based on finite differences (FDM) for simplicity. As shown in Tables 7 and 8, employing multiple levels achieves better performance under the same total number of smoothing steps, while also improving computational efficiency. Inference results are reported per sample under a relative tolerance of \(10^{-6}\).

6.7 Sensitivity to Points-Per-Wavelength↩︎

In this experiment, we evaluate the robustness of the proposed McMg architecture under different points-per-wavelength (\(\texttt{ppw}\)). The difficulty of solving the Helmholtz equation numerically is heavily dependent on the number of \(\texttt{ppw}\), we vary the wavenumber to sweep \(\texttt{ppw} \in \{6, 8, 10, 12, 14, 16\}\), transitioning from the challenging engineering limit to a well-resolved regime. The model is implemented using FDM discretization and trained on the OpenBreastUS breast dataset, following the protocols in Section 4.1.0.1.

The results, summarized in Table 9, reveal a clear correlation between grid resolution and convergence speed. In the well-resolved regime (\(\texttt{ppw} \ge 10\)), where the wave is smooth relative to the grid, the McMg converges rapidly, requiring only 4-5 iterations. As \(\texttt{ppw}\) decreases towards \(6\), the problem difficulty increases sharply due to the pollution effect and numerical dispersion inherent in standard discretizations. While classical geometric multigrid methods typically diverge in this coarse regime due to aliasing of oscillatory modes, McMg maintains stability, converging in 8.2 iterations at \(\texttt{ppw} \approx 6\). The increase in iterations at lower \(\texttt{ppw}\) indicates that the network must work harder to correct the solver’s coarse-grid dispersion mismatch, which sharpens as the wave becomes marginally resolved; this affects only the rate of convergence to the discrete solution, not the accuracy of that solution, and demonstrates the solver’s robustness on challenging, marginally resolved cases.

References↩︎

---

1. https://developer.nvidia.com/cudss↩︎