1 Introduction↩︎

Multiscale PDE problems, characterized by the coexistence of important features across multiple spatial and temporal scales, arise in a wide range of scientific and engineering domains. A representative example is given by Darcy Flow, \[\label{darcy95eq} \begin{align} - \nabla \cdot (\kappa \nabla u) &= f \quad &&\text{in } D, \\ u &= 0 \quad &&\text{on } \partial D, \end{align}\tag{1}\] where \(D \subset \mathbb{R}^d\), \(d = 2,3\) is a bounded physical domain, \(u : D \to \mathbb{R}\) denotes the pressure, \(f \in L^2(D)\) is a given source term and \(\kappa : D \to \mathbb{R}_{+}\) is a heterogeneous permeability coefficient. The Darcy Flow equation can be viewed as a heterogeneous extension of Poisson’s equation, which arises naturally across a range of physical settings involving transport through complex media. For example, it describes fluid flow in porous subsurface [1] formations, while analogous formulations govern heat conduction in materials [2] with spatially varying thermal properties and potential fields in electrically inhomogeneous media [3]. Across these problems, the spatial variability of the coefficient field plays a central role, as even fine-scale heterogeneities can significantly influence the global solution behavior, leading to pronounced multiscale effects. This interplay between fine-scale structure and macroscopic response poses significant challenges for both modeling and computation, often requiring high-resolution simulations to accurately capture the underlying physics. Such demands are particularly critical in practical engineering workflows, where repeated evaluations under varying material configurations or design parameters are often required. As a result, developing accurate and efficient solution strategies for this class of problems remains a key objective in both scientific computing and industrial applications [4], [5].

In this context, neural operator methods have emerged as promising alternatives, as they aim to learn mappings between infinite-dimensional function spaces and thereby approximate solution operators of parametric PDEs directly. However, despite their favorable computational properties and remarkable success in learning solution operators of parametric PDEs, neural operator models exhibit important limitations when applied to multiscale problems with rough and strongly heterogeneous inputs. Much of the existing neural operator literature focuses on relatively simplified flow configurations; for instance, benchmark studies involving Fourier Neural Operator (FNO) [6], Deep Operator Network (DeepONet) [7], and Transformer-based architectures such as Transolver [8], and related approaches frequently consider two-phase Darcy flow problems with moderately varying permeability fields (Figure 1 (a)). While these settings provide useful testbeds for model development and evaluation for such models, they do not fully reflect the complexity encountered in realistic subsurface flow applications. In particular, multi-phase Darcy flow in highly heterogeneous and highly oscillatory media (Figure 1 (b)-(e)) remains largely unexplored in the neural operator community, despite its relevance in groundwater hydrology, petroleum engineering, and carbon sequestration.

1.0.0.1 Contributions

In this work, we address these gaps and make the following contributions.

• We identify a key shortcoming of existing neural operator models when applied to multiscale problems with rough, high-contrast coefficients. Using Darcy flow as a prototypical example, we demonstrate the challenges posed by highly heterogeneous and multiscale permeability fields through a series of carefully designed datasets inspired by applications in industry that exhibit complex fine-scale features.

• We propose LOD-MSNO, an operator learning model specifically designed for multiscale PDE problems that follows a coefficient-learning paradigm. The approach adopts a hybrid strategy by using the finite element method-based Localized Orthogonal Decomposition (LOD) approximation as a strong numerical prior and employing neural networks to learn the associated multiscale basis functions and expansion coefficients.

• We provide a theoretical analysis showing that the proposed model converges to the true numerical LOD solution. As our main theorem, we derive a precise error bound with respect to the true analytical solution, separating the LOD discretization error from the generalization error of the neural surrogate.

• We demonstrate that our method can outperform established neural operator benchmarks on the proposed multiscale datasets. These results validate the feasibility of our idea and highlight its potential for reliable and efficient surrogate modeling of multiscale PDE problems.

2 Related work↩︎

2.0.0.1 Multiscale numerical methods

The numerical approximation of PDEs with highly heterogeneous and multiscale coefficients has been extensively studied in computational mathematics. Classical multiscale methods [9] and numerical homogenization [10] techniques aim to derive effective coarse-scale models that capture fine-scale effects implicitly. While successful in settings with clear scale separation, their performance degrades when coefficients exhibit strong local variations or lack periodicity. Localized approaches, most notably the Localized Orthogonal Decomposition (LOD) [11] method, address these limitations by enriching coarse finite element spaces with localized correctors that encode fine-scale information. LOD methods provide rigorous error estimates for elliptic [12] and parabolic [13] problems with rough and high-contrast coefficients but require the solution of many localized fine-scale problems, which can be computationally expensive in many-query settings.

2.0.0.2 Physics-informed machine learning

Physics-informed learning approaches such as Physics-Informed Neural Networks (PINNs) [14] and the closely related Deep Galerkin Method [15] embed the governing PDEs into the training objective through residual minimization. These methods offer mesh-free approximations and flexibility with respect to geometry and boundary conditions. However, they often suffer from optimization issues and have to be retrained for every input instance.

2.0.0.3 Neural operators

Neural operator methods aim to learn mappings between function spaces directly. DeepONet [16] represents operators as linear combinations of a branch and trunk network, with extensions such as Physics-informed DeepONet (PI-DeepONet) [17] or HyperDeepONet [18]. A complementary class of methods is based on spectral representations. The Fourier Neural Operator (FNO) [6] and its variants, including Convolutional Neural Operators (CNO) [19], U-shaped Neural Operators (UNO) [20] and the physics-informed PINO [21], have demonstrated strong performance on parametric PDE benchmarks. However, their reliance on global spectral representations can limit their ability to resolve localized fine-scale features in highly heterogeneous and multiscale settings. Other work has explored attention-based and graph-based operator learning [22]–[24]; for example, Transolver [8] is a representative example that leverages Transformer architectures to model long-range dependencies.

2.0.0.4 Hybrid numerical and machine learning approaches

Recent work has explored hybrid approaches that combine classical numerical methods, such as finite element or spectral discretizations, with neural networks in order to leverage the strengths of both paradigms [25]. In these methods, the solution is represented in a fixed basis, and the learning task is to infer the corresponding expansion coefficients. For example, Finite element operator network (FEONet) [26] integrates neural networks with the continuous Galerkin finite element method by directly learning solution coefficients, enabling physics-aware training without requiring labeled data. Related spectral approaches, such as SCLON [27] and ULGNet [28], learn coefficients in spectral or Galerkin bases, further demonstrating the effectiveness of combining classical numerical schemes with neural networks as a hybrid model for PDEs. Other hybrid AI-based approaches have also demonstrated substantial computational speedups and improved flexibility, with successful applications reported in areas such as solid mechanics [29], [30], fluid dynamics [31], climate science [32], and biomedical systems [33].

3 Preliminaries↩︎

3.1 Problem setting and Finite Element approximation↩︎

Let \(D \subset \mathbb{R}^d\), \(d \in \{2,3\}\), be a bounded polygonal (or polyhedral) Lipschitz domain with boundary \(\partial D\). We consider the Darcy flow problem with homogeneous Dirichlet boundary conditions: Given a permeability field \(\kappa\), find \(u \colon D \to \mathbb{R}\) such that \[\label{eq:darcy95prelim} \begin{align} - \nabla \cdot (\kappa \nabla u) &= f \quad &&\text{in } D, \\ u &= 0 \quad &&\text{on } \partial D, \end{align}\tag{2}\] where \(f \in L^2(D)\) and \(\kappa \colon D \to \mathbb{R}^{d \times d}\). Throughout this work, we assume that \(\kappa\) is uniformly elliptic and almost everywhere bounded. \[0 < \operatorname*{ess\,inf}_{x\in D} \kappa(x) \leq \kappa(x) \leq \operatorname*{ess\,sup}_{x\in D} \kappa(x) < \infty \qquad \text{for a.e. } x\in D .\] However, no further smoothness or scale separation assumptions on \(\kappa\) are imposed; in particular, \(\kappa\) may exhibit highly oscillatory or high-contrast behavior, i.e. the contrast of the permeability field \(\kappa\) defined by \[\operatorname{contrast}(\kappa) := \frac{ \operatorname*{ess\,sup}_{x\in D} \kappa(x) }{ \operatorname*{ess\,inf}_{x\in D} \kappa(x) }\] might have a very high value, which increases the difficulty of solving (2 ) numerically. For the approximation of this problem using FEM-based methods, we introduce two conforming simplicial triangulations of \(D\). Let \(\mathcal{T}_H\) denote a coarse triangulation with maximum element diameter \(H > 0\), and let \(\mathcal{T}_h\) be a fine triangulation with maximum element diameter \(h < H\). We assume that \(\mathcal{T}_h\) is a (possibly non-uniform) refinement of \(\mathcal{T}_H\) and that \(\overline{D} = \bigcup_{K \in \mathcal{T}_H} K = \bigcup_{K_h \in \mathcal{T}_h} K_h\). The fine mesh \(\mathcal{T}_h\) is chosen sufficiently small to resolve all relevant variations of the diffusion coefficient \(\kappa\) and is therefore referred to as a fine-scale discretization. We define the standard conforming \(P_1\) finite element spaces \[P_1(\mathcal{T}) := \left\{ v \in C^0(D) \;\middle|\; v|_K \text{ is affine for all } K \in \mathcal{T} \right\},\] and set \[V_h := P_1(\mathcal{T}_h) \cap H_0^1(D), \quad V_H := P_1(\mathcal{T}_H) \cap H_0^1(D),\] with \(V_H \subset V_h\). Let \(\mathcal{N}_H = \{ Z_i \mid 0 \le i \le N_H - 1 \}\) and \(\mathcal{N}_h = \{ z_j \mid 0 \le j \le N_h - 1 \}\) denote the sets of interior vertices (nodes) of the coarse and fine meshes, respectively, where \(N_H\) and \(N_h\) are their cardinalities. We denote by \[\{ \Phi_H^i \}_{i=1}^{N_H} \subset V_H, \qquad \{ \phi_h^j \}_{j=1}^{N_h} \subset V_h,\] the associated nodal basis functions. A fine-scale finite element approximation is obtained by seeking \(u_h \in V_h\) such that \[\label{eq:fine95fem} a(u_h, v_h) = \int_\Omega f v_h \, dx \quad \text{for all } v_h \in V_h.\tag{3}\] Writing \(u_h = \sum_{i=0}^{N_h-1} \boldsymbol{u}_i \phi_h^i\), this problem is equivalent to solving the linear system \(\mathbf{A}\boldsymbol{u} = \mathbf{f}\) with \(\boldsymbol{A}_{ij} = a(\phi_h^i,\phi_h^j)\) and \(\boldsymbol{f}_i = \int_\Omega f \phi_h^i , dx\). As already indicated, in order to accurately resolve the fine-scale features of a potentially highly-varying and oscillatory coefficient \(\kappa\), the mesh size \(h\) must typically be chosen very small, which may result in prohibitively high computational cost for solving the linear system. On the other hand, a coarse-scale finite element approximation based on \(V_H\) generally fails to capture the influence of the unresolved fine-scale variations in \(\kappa\) and therefore leads to poor accuracy.

3.2 Localized orthogonal decomposition (LOD)↩︎

The LOD method, originally proposed in [12], is a finite element method (FEM) designed to overcome this difficulty by constructing problem-adapted multiscale basis functions that combine the computational efficiency of coarse discretizations with the accuracy of fine-scale resolution.

3.2.1 Orthogonal decomposition of scales↩︎

The central idea of the LOD method is to construct an accurate approximation of the solution in a low-dimensional space of the same dimension as a coarse finite element space \(V_H\), while enriching this space with fine-scale information from a high-dimensional space \(V_h\). To this end, first define a projection (or quasi-interpolation) operator \(\mathcal{I}_H \colon V \to V_H\), which projects functions from the energy space \(V = H_0^1(D)\) onto its associated counterpart in the coarse finite element space \(V_H\). The operator \(\mathcal{I}_H\) is assumed to be a projection in the algebraic sense, i.e., \[\mathcal{I}_H \circ \mathcal{I}_H = \mathcal{I}_H,\] and to satisfy appropriate stability properties [11]. Several realizations of such a projection are possible. Typical examples include the global \(L^2\)-orthogonal projection onto the coarse finite element space \(V_{H}\), defined by \[(\mathcal{I}_H v, \Phi_H)_{L^2(\Omega)} = (v, \Phi_H)_{L^2(\Omega)} \quad \text{for all } \Phi_H \in V_{H},\] as well as locally defined \(L^2\)-projections based on nodal patches. These local constructions exploit projections onto polynomial spaces restricted to coarse neighborhoods and lead to operators that satisfy the required stability properties on quasi-uniform meshes. In contrast, nodal interpolation operators generally fail to provide sufficient stability in the \(H^1\)-norm and are therefore not suitable for the LOD framework; see [11] for a detailed discussion of admissible choices. We define the space of microscopic fine-scale functions as \[W := \ker(\mathcal{I}_H) = \left\{ v \in V \;\middle|\; \mathcal{I}_H(v) = 0 \right\}.\] The space \(W\) contains all functions that encode fine-scale features that cannot be represented in the coarse space \(V_H\) and yields an orthogonal decomposition \(V = V_H \oplus W\) w.r.t. the standard inner product on \(V\). In order to construct an improved approximation space compared to \(V_H\), we exploit the problem-specific inner product on \(V\) induced by the bilinear form \[a(u,v) := \int_D \kappa \nabla u \cdot \nabla v \, dx.\] and define the improved multiscale approximation space as \[V_H^{\mathrm{ms}} := \left\{ v_H^{\mathrm{ms}} \in V \;\middle|\; a(v_H^{\mathrm{ms}}, w) = 0 \;\; \text{for all } w \in W \right\},\] which is nothing else but the orthogonal complement of the subspace \(W\) in \(V\) with respect to the inner product \(a\), leading to the \(\mathcal{A}\)-orthogonal decomposition \(V = V_H^{\mathrm{ms}} \oplus W\). By construction, \(V_H^{\mathrm{ms}}\) has the same dimension as \(V_H\) and consists of functions that are \(\mathcal{A}\)-orthogonal to all unresolved fine-scale components. Based on this decomposition, we define the ideal numerical homogenization method as follows: given \(f \in L^2(D)\), find \(u_H^{\mathrm{ms}} \in V_H^{\mathrm{ms}}\) such that \[a(u_H^{\mathrm{ms}}, v) = \int_D f v \, dx \quad \text{for all } v \in V_H^{\mathrm{ms}}.\] It is ideal in the sense that \(V^{ms}_H\) contains all coarse-resolvable fine-scale features of the multiscale coefficient \(\kappa\).

3.2.2 Practical realization via localization↩︎

3.2.2.1 Finite-dimensional realization of the LOD method

The ideal numerical homogenization method introduced above is formulated in the infinite-dimensional multiscale space \(V_H^{\mathrm{ms}} \subset V\). In order to make this formulation computationally feasible, we next clarify how the corresponding variational problem can be solved in practice. To this end, consider the correction operator \(Q \colon V \to W\), defined as the \(a(\cdot,\cdot)\)-orthogonal projection onto the detail space \(W = \ker(\mathcal{I}_H)\). More precisely, for a given \(v \in V\), the correction \(Qv \in W\) is defined as the unique solution of \[\label{correction95Def} a(Qv, w) = a(v, w) \quad \text{for all } w \in W.\tag{4}\] By construction, \((1 - Q)v\) is \(a\)-orthogonal to \(W\). One can show that the operator \((1 - Q) \colon V_H \to V_H^{\mathrm{ms}}\) is bijective, with inverse given by the projection \(\mathcal{I}_H\). Consequently, every function \(v_H^{\mathrm{ms}} \in V_H^{\mathrm{ms}}\) can be written uniquely as \(v_H^{\mathrm{ms}} = (1 - Q) v_H\) for some \(v_H \in V_H\). Using this observation, the ideal multiscale problem can be reformulated as follows: find \(u_H \in V_H\) such that \[a\bigl((1 - Q) u_H, (1 - Q) v_H\bigr) = \int_D f (1 - Q) v_H \, dx\] for all \(v_H \in V_H\). This formulation is more favorable from a computational perspective, as it depends on the discrete, low-dimensional unknown \(u_H \in V_H\) instead of the infinite-dimensional object \(u_H^{\mathrm{ms}} \in V_H^{\mathrm{ms}}\).

3.2.2.2 Practical computation via localization

In practice, the numerical construction of the correction operator \(Q\) itself is localized to small patches around each coarse element in order to make the method computationally feasible. This localization is motivated by the exponential decay of the corrector functions away from their associated coarse elements; see Chapter 4 in [11]. For a given coarse element \(T \in \mathcal{T}_H\), we define a patch \(\mathcal{N}^\ell(T)\) consisting of all coarse elements within \(\ell\) layers of \(T\) recursively by \[\begin{align} \mathcal{N}^0(T) &:= T, \\ \mathcal{N}^\ell(T) &:= \bigcup \left\{ T_H \in \mathcal{T}_H \;\middle|\; T_H \cap \mathcal{N}^{\ell-1}(T) \neq \emptyset \right\}, \quad \ell \in \mathbb{N}. \end{align}\] and compute the localized correction operator \(Q_h^\ell\) by solving fine-scale elliptic problems of the form 4 on these patches. To be precise, let \(W^\ell(T) \subset W\) be the space consisting of functions in \(W\) supported on \(\mathcal{N}^\ell(T)\) and extended by zero outside this region. The localized element corrector \[Q_T^\ell : V \rightarrow W^\ell(T)\] is then defined as the solution of \[a(Q_T^\ell v, w) = \int_T \kappa \nabla v \cdot \nabla w \, dx \quad \text{for all } w \in W^\ell(T). \label{eq:localized95element95corrector}\tag{5}\]

For practical computations, all spaces are discretized using the fine-scale finite element space \(V_h\). For a given \(v_h \in V_h\), the fully discrete localized element corrector \[Q_{T,h}^\ell(v_h) \in W^\ell(T) \cap V_h\] is defined as the solution of \[\label{eq:discrete95localized95corrector} a\bigl(Q_{T,h}^\ell(v_h), w_h\bigr) = \int_T \kappa \nabla v_h \cdot \nabla w_h \, dx \quad \text{for all } w_h \in W^\ell(T) \cap V_h.\tag{6}\] The corresponding discrete localized correction operator is then given by \[Q_h^\ell := \sum_{T \in \mathcal{T}_H} Q_{T,h}^\ell.\]

Finally, the fully discrete Localized Orthogonal Decomposition method reads as follows: given \(f \in L^2(D)\), find \(u_{H,\ell,h}^{\mathrm{ms}} \in V_{H,\ell}^{\mathrm{ms},h}\) such that \[\label{eq:variational95lod95method} a(u_{H,\ell,h}^{\mathrm{ms}}, v) = \int_D f v \, dx \quad \text{for all } v \in V_{H,\ell}^{\mathrm{ms},h},\tag{7}\] where the multiscale space is defined by \[V_{H,\ell}^{\mathrm{ms},h} := \operatorname{span} \left\{ \Phi_H^i - Q_h^\ell(\Phi_H^i) \;\middle|\; Z_i \in \mathcal{N}_H \right\}.\] This space has the same dimension as the coarse finite element space \(V_H\), and the method can be interpreted as a Galerkin approximation of the original problem in a problem-adapted multiscale space. Both the correction operator \(Q_h^\ell\) and the final solution \(u_{H,\ell,h}^{\mathrm{ms}}\) can be obtained by assembling and solving the corresponding linear systems; see [34] for implementation details. In what follows, we shall omit the coarse and fine mesh sizes \(H\) and \(h\) and the localization parameter \(\ell\) and call the solution \(u_{H,\ell,h}^{\mathrm{ms}}\) to (7 ) simply \(u_{\mathrm{LOD}}\).

3.3 Limitations of the LOD method↩︎

Despite its effectiveness for multiscale problems without any assumption about the periodicity of \(\kappa\), the LOD method suffers from notable computational limitations. In particular, LOD-based discretization requires the solution of \(|\mathcal{T}_H|\) fine-scale auxiliary problems in overlapping subdomains \(\mathcal{N}^\ell(T)\) for each realization of the random coefficient. As the localization parameter \(\ell\) increases, these subdomains grow in size and overlap more strongly, leading to a rapidly increasing computational cost for the construction of the correction operator \(Q_h^{\ell}\). However, theoretical accuracy results typically require choosing \(\ell \approx C \lvert \log H \rvert\) for some constant \(C>0\), which may result in relatively large localization patches in practice. This renders the computation of the correctors a major computational bottleneck.

Moreover, the corrected basis functions lose the local support of standard finite element nodal basis functions: Instead of being supported on a single coarse element \(T\), the LOD basis functions have support on the entire patch \(\mathcal{N}^\ell(S)\), thereby significantly polluting the sparsity structure of the resulting system matrix. This effect is further exacerbated in settings where the resulting sparsity-polluted linear systems must be solved repeatedly. These considerations motivate the development of fully surrogate-based models that bypass both the explicit computation of the basis corrections \(Q_{h}^{\ell}\) and the solution of the associated dense or weakly sparse linear systems.

4 Proposed method: LOD Multiscale Neural Operator (LOD-MSNO)↩︎

The overall goal of LOD-MSNO is to learn the solution operator \[\mathcal{G}: \kappa \mapsto u,\] for the multiscale Darcy problem by using neural networks to surrogate the computationally expensive components of the LOD method. We assume the permeability is a random field \(\kappa=\kappa(\omega)\) on a probability space \((\Omega,\mathcal{F},\mathbb{P})\), where \(\Omega\) is assumed to be compact and \(\omega\in\Omega\) parameterizes one realization of the coefficient field. In the experiments, \(\kappa\) may be generated from Gaussian or lognormal random fields; the specific input distributions we used are described in Section 5. Accordingly, both the input permeability field and solution are interpreted as Bochner-type functions written as \(\kappa(x,\omega)\) and \(u(x,\omega)\) for \((x,\omega) \in D \times \Omega\).

Central to the LOD method are the corrected FEM basis functions \[\label{eq:lod95basis} \Psi_H^j(\omega) := \Phi_H^j - Q_{h,\kappa(\omega)}^{\ell}(\Phi_H^j), \qquad j = 1,\dots,N_H,\tag{8}\] where \(\{\Phi_H^j\}\) denotes the standard coarse nodal basis functions and \(Q_{h,\kappa(\omega)}^{\ell}\) is the localized correction operator computed from fine-scale problems posed on overlapping patches of diameter proportional to the localization parameter \(\ell\). The dependence on \(\omega\) enters through the bilinear form \[a_{\omega}(v,w) := \int_D \kappa(x,\omega)\nabla v\cdot\nabla w\,dx,\] which defines the localized corrector problems. Consequently, the corrected basis functions \(\Psi_H^j(\omega)\) are problem-adapted and depend on the realization of the coefficient field. Using the corrected basis \(\{\Psi_H^j(\omega)\}\), the LOD solution is represented as \[\label{lod95linearcombi} u_{\mathrm{LOD}}(x,\omega) = \sum_{j=1}^{N_H} \boldsymbol{u}_{\mathrm{LOD}}^{(j)}(\omega) \Psi_H^j(x,\omega),\tag{9}\] where the coefficient vector \(\boldsymbol{u}_{\mathrm{LOD}}(\omega) \in \mathbb{R}^{N_H}\) is obtained by solving the linear system \[\label{eq:lod95linear95eq} \boldsymbol{A}_{\mathrm{LOD}}(\omega) \boldsymbol{u}_{\mathrm{LOD}}(\omega) = \boldsymbol{f}_{\mathrm{LOD}}(\omega).\tag{10}\] Here, \(\boldsymbol{A}_{\mathrm{LOD}}(\omega)\) is the LOD stiffness matrix and \(\boldsymbol{f}_{\mathrm{LOD}}(\omega)\) is the LOD load vector, with entries \[{\boldsymbol{A}_{\mathrm{LOD}}(\omega)}_{ij} := a_{\omega}(\Psi_H^i(\omega), \Psi_H^j(\omega)), \qquad {\boldsymbol{f}_{\mathrm{LOD}}(\omega)}_i := \int_D f\Psi_H^i(x,\omega)\,dx.\] LOD-MSNO targets this parameter-to-solution map by predicting the coefficient vector \(\boldsymbol{u}_{\mathrm{LOD}}(\omega)\) directly with a neural network, while avoiding explicit assembly and solution of 10 at inference time. At the same time, surrogating the construction of the correction operator \(Q_{h,\kappa(\omega)}^{\ell}\)—and thereby obtaining the corrected multiscale basis functions \(\Psi_H^i(\omega)\)—poses a separate and equally important challenge.

4.1 LOD-MSNO (Coeff): Learning the coefficients↩︎

The first component of our surrogate focuses on learning the coefficient map \[\omega \mapsto \boldsymbol{u}_{\mathrm{LOD}}(\omega),\] or equivalently \(\kappa(\cdot,\omega)\mapsto \boldsymbol{u}_{\mathrm{LOD}}(\omega)\), without assembling and solving 10 explicitly at inference time. To train this neural network surrogate, we draw i.i.d. realizations \(\omega_i\in\Omega\) and evaluate the corresponding permeability fields \(\kappa_i:=\kappa(\cdot,\omega_i)\). We then construct a dataset of tuples \[\{( \boldsymbol{A}_{\mathrm{LOD}}(\omega_i), \boldsymbol{u}_{\mathrm{LOD}}(\omega_i) )\}_{i=1}^M,\] where \(\boldsymbol{A}_{\mathrm{LOD}}(\omega_i)\) denotes the LOD stiffness matrix and \(\boldsymbol{u}_{\mathrm{LOD}}(\omega_i)\) denotes the ground-truth coefficient vector for the realization \(\kappa(\cdot,\omega_i)\). Since the space \(V_{H,\ell}^{\mathrm{ms},h}(\omega_i) = \operatorname{span} \left\{ \Psi_H^j(\omega_i) \right\}\) is typically low-dimensional, storing the matrices \(\boldsymbol{A}_{\mathrm{LOD}}(\omega_i)\) remains feasible in terms of memory consumption. Moreover, these matrices have to be computed during data generation to obtain the ground-truth coefficients \(\boldsymbol{u}_{\mathrm{LOD}}(\omega_i)\). Rather than employing a standard \(L^2\)- or mean-square-based regression loss, we adopt an energy-based loss that reflects the variational structure of the underlying elliptic problem. At the population level, for a predicted coefficient map \(\boldsymbol{u}:\Omega\to\mathbb{R}^{N_H}\), we define \[\label{eq:energy95loss95population} \mathcal{L}_{E}(\boldsymbol{u}) := \left\| \boldsymbol{A}_{\mathrm{LOD}}^{1/2}(\omega) \bigl(\boldsymbol{u}(\omega)-\boldsymbol{u}_{\mathrm{LOD}}(\omega)\bigr) \right\|_{L^1(\Omega)}.\tag{11}\] In practice, this expectation over the probability space is replaced by a Monte Carlo approximation using i.i.d. samples \(\omega_i\sim\mathbb{P}\), which yields the empirical loss \[\label{eq:energy95loss} \mathcal{L}_{E,M}(\boldsymbol{u}) := \frac{1}{M}\sum_{i=1}^M \left| \boldsymbol{A}_{\mathrm{LOD}}^{1/2}(\omega_i) \bigl(\boldsymbol{u}(\omega_i)-\boldsymbol{u}_{\mathrm{LOD}}(\omega_i)\bigr) \right|.\tag{12}\]

The choice of the energy-based loss 12 is directly motivated by the variational formulation of the underlying elliptic PDE. At the discrete level, this variational principle is associated with the (discrete) energy norm of a coefficient vector \(\boldsymbol{v}\), which is given by \[\|\boldsymbol{v}\|_{\boldsymbol{A}_{\mathrm{LOD}}}^2 = \boldsymbol{v}^{\top} \boldsymbol{A}_{\mathrm{LOD}} \boldsymbol{v}.\] Consequently, the proposed loss 12 measures the error between the predicted and true coefficients in the natural energy norm associated with the elliptic operator. We hypothesize that employing this energy-based loss can be beneficial for generalization, as it is more closely aligned with the variational principle underlying the PDE. We empirically observe that training with the energy-based loss 12 leads to improved accuracy and better generalization, particularly for coefficients exhibiting strong heterogeneity and high contrast. An ablation study comparing different loss functions is provided in Section 5.2.2. In the next section, we develop an error estimate for the LOD-MSNO (Coeff) model and analyze how using the energy-based loss affects its convergence properties.

4.2 Error analysis for LOD-MSNO (Coeff)↩︎

Throughout this subsection, we assume that the corrected LOD basis functions \(\{\Psi_H^j(\omega)\}_{j=1}^{N_H}\) are available for each realization \(\omega\in\Omega\). Hence, the basis functions are treated as given, and the learning task reduces to approximating the coefficient map associated with the LOD representation. Let \(\mathcal{N}\) be any class of \(\mathbb{R}^{N_H}\)-valued functions parameterized by neural networks. We distinguish the following objects:

• \(u(x,\omega)\) denotes the true solution of the Darcy problem 1 .

• \(u_{\mathrm{LOD}}(x,\omega)\) denotes the numerical LOD solution obtained as a linear combination of the LOD basis by solving the linear system 10 .

• \(\widehat{u}_{\mathrm{LOD}}^{E}(x,\omega)\) denotes the solution reconstructed from coefficients predicted by a network trained with the population loss \(\mathcal{L}_{E}\) in 11 .

• \(\widehat{u}_{\mathrm{LOD}}^{E,M}(x,\omega)\) denotes the actual LOD-MSNO (Coeff) prediction reconstructed from coefficients trained with the empirical Monte Carlo loss \(\mathcal{L}_{E,M}\) in 12 .

More precisely, if \[\widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E} \in \operatorname*{argmin}_{\boldsymbol{u}\in\mathcal{N}}\mathcal{L}_{E}(\boldsymbol{u}), \qquad \widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,M} \in \operatorname*{argmin}_{\boldsymbol{u}\in\mathcal{N}}\mathcal{L}_{E,M}(\boldsymbol{u}),\] then the corresponding reconstructed LOD solutions are \[\begin{align} \widehat{u}_{\mathrm{LOD}}^{E}(x,\omega) &:= \sum_{j=1}^{N_H} \widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,(j)}(\omega)\Psi_H^j(x,\omega),\\ \widehat{u}_{\mathrm{LOD}}^{E,M}(x,\omega) &:= \sum_{j=1}^{N_H} \widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,M,(j)}(\omega)\Psi_H^j(x,\omega). \end{align}\] In particular, throughout this subsection, bold symbols denote coefficient vectors, whereas non-bold symbols denote the corresponding (reconstructed) functions. The goal of the following analysis is to quantify the discrepancy between the true solution and the learned LOD-MSNO (Coeff) prediction. To this end, we decompose the total error into three contributions: \[\label{eq:error95decomposition} \begin{align} u(x,\omega)-\widehat{u}_{\mathrm{LOD}}^{E,M}(x,\omega) ={}& \underbrace{u(x,\omega)-u_{\mathrm{LOD}}(x,\omega)}_{\text{(I) LOD approximation error}} \\ &+\underbrace{u_{\mathrm{LOD}}(x,\omega)-\widehat{u}_{\mathrm{LOD}}^{E}(x,\omega)}_{\text{(II) population learning error}} \\ &+\underbrace{\widehat{u}_{\mathrm{LOD}}^{E}(x,\omega)-\widehat{u}_{\mathrm{LOD}}^{E,M}(x,\omega)}_{\text{(III) empirical/generalization error}}. \end{align}\tag{13}\] The three terms correspond to the LOD discretization error, the approximation error induced by the neural network class, and the generalization error arising from finite sampling. The first term is controlled by the classical LOD approximation theory. We state the corresponding result here:

The proof of the above result can be found in [11]. In particular, term (I) in 13 converges like \(H^2\) in the \(L^2(D)\) norm when \(\ell\gtrsim |\log H|\). In the following, we now develop bounds for terms (II) and (III) in Theorems 2 and 3, respectively. Their proofs are given in the appendix. We first start with a bound for (II).

Next, we shall give a bound for term (III).

If the neural network class \(\mathcal{N}\) is chosen appropriately through the architecture, such that the associated Rademacher complexity satisfies \(R_M(\mathcal{G}^{E})\to0\) as \(M\to\infty\), and if \(\mathcal{N}\) is sufficiently expressive to approximate the LOD coefficient map, then the learned coefficient surrogate converges to the ground-truth LOD coefficients \(\boldsymbol{u}_{\mathrm{LOD}}\). Since the final LOD-MSNO(Coeff) prediction is obtained by reconstructing the physical solution from these coefficients and the multiscale basis functions, it remains to transfer this coefficient convergence to convergence of the reconstructed solution. The following theorem establishes this result and is our main theoretical result. It provides a precise error bound between the exact solution of equation 1 and the LOD-MSNO (Coeff) coefficient surrogate, combining the LOD discretization error with the approximation and generalization errors of the neural coefficient model.

Theorem 4 establishes that the total error between the true solution of (1 ) is primarily governed by three contributions. The first term reflects the approximation properties of the LOD method and recovers the optimal bound when \(\ell \approx |\log(H)|\). The second term captures the generalization capability of the neural network class, which can be quantified via its Rademacher complexity. The third term depends on both the expressive power of the neural network class \(\mathcal{N}\) and the spectral properties of the LOD stiffness matrix corresponding to a given input \(\omega \in \Omega\). By picking a neural network class \(\mathcal{N}\) with strong enough approximation (third term) and generalization properties (second term), and by choosing \(H\) small enough (first term), we can achieve that the expected error between the true solution \(u\) and the LOD-MSNO (Coeff) prediction \(\widehat{u}_{\mathrm{LOD}}^{E,M}\) converges to zero.

4.3 LOD-MSNO (Joint): Learning the basis functions and the combined model↩︎

The coefficient model described so far assumes that the corrected LOD basis functions \(\{\Psi_H^j(\omega)\}_{j=1}^{N_H}\) are available for every input realization. However, computing these basis functions by solving localized fine-scale corrector problems of the form (6 ) for every patch constitutes the major computational bottleneck of the classical LOD method. In order to remove this remaining bottleneck, we aim to also learn the coefficient-dependent basis map \[\label{eq:basis95map95learning} \kappa(\cdot,\omega) \longmapsto \{\Psi_H^j(\cdot,\omega)\}_{j=1}^{N_H}.\tag{14}\] that maps an input permeability field \(\kappa\) to its corresponding set of LOD-corrected basis functions. Rather than learning the multiscale basis functions directly, we exploit their LOD construction. For each coarse basis function \(\Phi_H^j\), the corresponding localized multiscale basis function is given by \[\label{eq:learned95basis95via95corrector} \Psi_H^j(\omega) := \Phi_H^j - Q_h^\ell(\Phi_H^j)(\omega) = \Phi_H^j - \sum_{T\in\mathcal{T}_H} Q_{T,h}^\ell(\Phi_H^j)(\omega).\tag{15}\] Consequently, it is sufficient to learn the local correction operators \(Q_{T,h}^\ell\) on each patch \(\mathcal{N}^\ell(T)\).

We now rewrite the discrete localized corrector problem in a form that is suitable for supervised learning. Let \(T\in\mathcal{T}_H\) be fixed and let \(\{\phi_{\ell,T,m}\}_{m=1}^{N_{\ell,h}(T)}\) denote the fine-scale nodal basis functions associated with the fine interior degrees of freedom in the patch \(\mathcal{N}^\ell(T)\). For a given coarse basis function \(\Phi_H^j\), we write the local corrector as \[\label{eq:local95corrector95coeff95expansion} Q_{T,h}^\ell(\Phi_H^j)(x,\omega) = \sum_{m=1}^{N_{\ell,h}(T)} w_{T,j,m}(\omega)\,\phi_{\ell,T,m}(x),\tag{16}\] and collect the local correction coefficients in the vector \(\boldsymbol{w}_{T,j}(\omega)\in\mathbb{R}^{N_{\ell,h}(T)}\). The defining local corrector equation is \[\label{eq:discrete95localized95corrector95learning} a_{\omega}\bigl(Q_{T,h}^\ell(\Phi_H^j), w_h\bigr) = \int_T \kappa(x,\omega)\nabla \Phi_H^j \cdot \nabla w_h \, dx \quad \text{for all } w_h \in W^\ell(T) \cap V_h.\tag{17}\] The constraint \(Q_{T,h}^\ell(\Phi_H^j)\in W^\ell(T)=\ker(\mathcal{I}_H)\cap \{v:\operatorname{supp}(v)\subset\mathcal{N}^\ell(T)\}\) can be imposed algebraically by a local constraint matrix \(\boldsymbol{C}_{T}\) representing the action of \(\mathcal{I}_H\) on the patch degrees of freedom. More precisely, if \(\boldsymbol{A}_{T}(\omega)\in\mathbb{R}^{N_{\ell,h}(T)\times N_{\ell,h}(T)}\) and \(\boldsymbol{r}_{T,j}(\omega)\in\mathbb{R}^{N_{\ell,h}(T)}\) are defined by \[\begin{align} \label{eq:local95corrector95matrix95entries} (\boldsymbol{A}_{T}(\omega))_{mn} &:= \int_{\mathcal{N}^\ell(T)} \kappa(x,\omega)\nabla\phi_{\ell,T,n}(x)\cdot \nabla\phi_{\ell,T,m}(x)\,dx, \\ (\boldsymbol{r}_{T,j}(\omega))_{m} &:= \int_T \kappa(x,\omega)\nabla\Phi_H^j(x)\cdot \nabla\phi_{\ell,T,m}(x)\,dx, \end{align}\tag{18}\] then the local correction coefficients are obtained from the saddle-point system \[\label{eq:local95corrector95saddle95point} \begin{pmatrix} \boldsymbol{A}_{T}(\omega) & \boldsymbol{C}_{T}^{\top} \\ \boldsymbol{C}_{T} & \boldsymbol{0} \end{pmatrix} \begin{pmatrix} \boldsymbol{w}_{T,j}(\omega) \\ \boldsymbol{\lambda}_{T,j}(\omega) \end{pmatrix} = \begin{pmatrix} \boldsymbol{r}_{T,j}(\omega) \\ \boldsymbol{0} \end{pmatrix}.\tag{19}\]

Here, \(\boldsymbol{\lambda}_{T,j}(\omega)\) is a vector of Lagrange multipliers that enforces the fine-scale constraint \(\mathcal{I}_H Q_{T,h}^\ell(\Phi_H^j)=0\) on the patch. If we let \(\boldsymbol{w}_T \in \mathbb{R}^{c_d \times N_{\ell,h}}\) be correction values associated with the patch centered around \(T \in \mathcal{T}_H\) for each of the \(c_d\) degrees of freedom in the element \(T\), the numerical construction of the LOD basis produces, for each realization \(\omega\), a collection of local coefficient matrices \[\label{eq:local95corrector95tensor} \boldsymbol{W}(\omega) := \bigl\{\boldsymbol{w}_{T}(\omega) \;:\; T\in\mathcal{T}_H\}.\tag{20}\] Equivalently, after choosing a fixed ordering of patches and local fine degrees of freedom, this collection can be viewed as a sparse tensor whose entries are the numbers \(w_{T,j,m}(\omega)\) in 16 . Based on the above formulation, instead of trying to directly approximate (14 ), our basis-correction network is trained to approximate the map \[\label{eq:correction95map95learning} \kappa(\cdot,\omega) \longmapsto \widehat{\boldsymbol{W}}(\omega)\tag{21}\]

Given a training set \[\label{eq:basis95training95dataset} \mathcal{D}_{B} := \bigl\{(\kappa_i,\boldsymbol{W}_i)\bigr\}_{i=1}^{M}, \qquad \boldsymbol{W}_i=\boldsymbol{W}(\omega_i),\tag{22}\] where \(\boldsymbol{W}_i\) is computed by solving 19 for all relevant patches and coarse basis functions, we use the supervised loss function \[\label{eq:basis95supervised95loss} \mathcal{L}_{basis,M}(\theta_B) := \frac{1}{M} \sum_{i=1}^{M} \sum_{T\in\mathcal{T}_H} \sum_{j=1}^{N_H} \frac{1}{N_{\ell,h}(T)} \left\| \boldsymbol{P}_T\widehat{\boldsymbol{w}}_{T,j}(\omega_i) - \boldsymbol{w}_{T,j}(\omega_i) \right\|.\tag{23}\] the learn the parameters \(\theta_B\) of the basis-correction network. Here, \(\boldsymbol{P}_T\) denotes the nullspace projector \[\boldsymbol{P}_T = I - \boldsymbol{C}_T^T {(\boldsymbol{C}_T\boldsymbol{C}_T^T)}^{-1}\boldsymbol{C}_T\] that ensures the predicted coefficients \(\widehat{\boldsymbol{w}}_{T,j}\) are in the space \(\ker(\mathcal{I}_H)\cap \{v:\operatorname{supp}(v)\subset\mathcal{N}^\ell(T)\}\). In practice, the summation over \(j\) only involves coarse basis functions whose support intersects the element \(T\), since all other right-hand sides \(\boldsymbol{r}_{T,j}(\omega)\) vanish. Once the local corrections have been predicted, the learned multiscale basis is assembled by \[\label{eq:learned95multiscale95basis} \widehat{\Psi}_{H}^j(\omega) := \Phi_H^j - \sum_{T\in\mathcal{T}_H} \sum_{m=1}^{N_{\ell,h}(T)} \widehat{\boldsymbol{w}}_{T,j}^{(m)}(\omega)\,\phi_{\ell,T,m}.\tag{24}\] This gives a learned approximation of the full LOD basis map in 14 without solving local saddle-point systems at inference time. To obtain a single neural operator from the permeability field to the fine-scale solution, both components are combined through the LOD reconstruction formula. Let \(\widehat{\boldsymbol{u}}_{\mathrm{LOD}}(\omega)\) denote the coefficient prediction and let \(\widehat{\boldsymbol{W}}(\omega)\) denote the predicted local correction coefficients. LOD-MSNO (Joint) reconstructs the full solution as \[\label{eq:full95solution95pred} \begin{align} \widehat{u}_{LOD-MSNO} &= \sum_{j=1}^{N_H} \widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,M,(j)}(\omega)\widehat{\Psi}_H^j(x,\omega) \\ &= \sum_{j=1}^{N_H} \widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,M,(j)}(\omega) \left( \Phi_H^j - \sum_{T\in\mathcal{T}_H} \widehat{Q}_{T,h}^\ell(\Phi_H^j)(\omega) \right) \\ & = \sum_{j=1}^{N_H} \widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,M,(j)}(\omega)\Phi_H^j - \sum_{j=1}^{N_H}\widehat{\boldsymbol{u}}_{\mathrm{LOD}}^{E,M,(j)}(\omega)\sum_{T\in\mathcal{T}_H} \sum_{m=1}^{N_{\ell,h}(T)} \widehat{\boldsymbol{w}}_{T,j}^{(m)}(\omega)\,\phi_{\ell,T,m} \end{align}\tag{25}\]

The coefficient network and the basis-correction network define two complementary surrogates. The first network approximates the coefficient map \(\kappa(\cdot,\omega)\mapsto\boldsymbol{u}_{\mathrm{LOD}}(\omega)\), whereas the second network approximates the correction map \(\kappa(\cdot,\omega)\mapsto\boldsymbol{W}(\omega)\). LOD-MSNO then predicts the full solution via (25 ). The full procedure to train the full LOD-Multiscale Neural Operator is summarized in Algorithm 3.

5 Experiments↩︎

We evaluate the proposed LOD-MSNO on the steady Darcy flow (1 ) problem, where the objective is to approximate the LOD representation of the solution from the permeability coefficient \(\kappa\). As already indicated through naming, we distinguish two models. LOD-MSNO (Coeff) is the coefficient model that assumes the exact corrected multiscale basis \(\{\Psi_H^i\}\) to be available and learns the map from \(\kappa\) to the expansion coefficients \(\boldsymbol{u}_{\mathrm{LOD}}\) in 9 . The full PDE solution is then reconstructed as a linear combination of the exact corrected basis functions with the learned coefficients. LOD-MSNO (Joint) is the combination of the coefficient model and the basis model that predicts, by contrast, both the coefficients \(\widehat{\boldsymbol{u}}_{\mathrm{LOD}}\) and the localized correction operators that define the coefficient-dependent multiscale basis as described in Section 4.3 and reconstructs the full solution via (25 ). Thus, the coefficient model only replaces the LOD-system (10 ) assembly and solve, whereas the basis model additionally replaces the localized basis corrector computations and predicts the solution to (1 ) without the need of any numerical computations.

5.0.0.1 Baseline models

We compare our method against a range of established baselines to ensure a balanced and meaningful evaluation. First, we include FNO [6], CNO [19], and UNO [20], which are widely used and representative neural operator architectures based on the iterative application of kernel integrals. Similarly, Transolver [8] is incorporated as a multi-purpose transformer-based operator model capable of capturing long-range dependencies. These models serve as standard data-driven operator learning benchmarks. To ensure fairness with respect to physics-informed approaches, we additionally compare against PINO [21], as our method also incorporates physical structure—specifically through the use of problem-adapted basis functions and the inclusion of the LOD stiffness matrix \(\boldsymbol{A}_{LOD}\) in the loss function. We further include FEONet [26] as a representative hybrid neural network–numerical method, since it also follows a coefficient-learning paradigm in which the solution is reconstructed as a linear combination of basis functions. In contrast to LOD-MSNO, however, FEONet uses fixed coarse \(P1\) finite element basis functions rather than coefficient-dependent corrected LOD basis functions. Finally, we consider PI-DeepONet [17], which is both physics-informed and structurally similar to our approach in that it represents the solution as a linear expansion over a learned set of basis functions. Our LOD-MSNO models are trained to predict the LOD coefficient vector \(\boldsymbol{u}_{LOD}\), from which the full fine-scale solution is reconstructed using the corrected LOD basis.

5.1 Benchmark datasets↩︎

As depicted in Figure 1, training and test datasets are generated using four distinct classes of permeability fields designed to probe different multiscale and high-contrast regimes. In the following, we refer to them as lognormal1, lognormal2, quantile, and checkerboard. These datasets are specifically constructed to highlight a key shortcoming of neural operator models when applied to multiscale problems and are inspired by industry-relevant real-world applications. A more detailed description of the dataset characteristics and its generation process can be found in 9.

5.1.0.1 Lognormal coefficients

For the lognormal datasets, we consider permeability fields of the form \(\kappa(x) = \exp(Z(x))\) where \(Z\) is a centered Gaussian random field with Whittle–Matérn covariance with smoothness parameter \(\nu>0\), correlation length \(\ell>0\), and variance \(\sigma^2>0\). We distinguish two regimes: (1) Lognormal1 with \(\sigma=1.0\), \(\ell=0.1\), (2) Lognormal2 with \(\sigma=2.0\), \(\ell=0.3\). The second setting exhibits significantly stronger oscillations and substantially higher contrast, making it particularly challenging for operator learning methods.

5.1.0.2 Quantile coefficients

For the quantile dataset, we first sample a Gaussian random field \(Z\) with squared-exponential covariance and define a lognormal field \(\exp(Z)\). The resulting values are then discretized into a finite set of prescribed permeability levels via empirical quantiles and randomly permuted across the domain, yielding piecewise-constant coefficients with sharp jumps and high contrast.

5.1.0.3 Checkerboard coefficients.

The checkerboard coefficients are defined as piecewise-constant random fields on a fine Cartesian grid, where each cell independently takes values from a fixed set of permeability levels. This yields coefficients with discontinuities at the grid scale and no spatial correlation.

5.2 Evaluation of LOD-MSNO↩︎

5.2.1 Comparison of baseline models and LOD-MSNO↩︎

In this section, we fix the coarse and fine mesh sizes to \(H=2^{-4}\) and \(h=2^{-7}\) for the evaluation of LOD-MSNO. LOD-MSNO (Coeff) uses the exact corrected multiscale basis functions and only replaces the LOD coefficient solve with a learned coefficient predictor. The coefficient model is realized by GlobalCoeffNet, which maps the normalized permeability field to the coefficient vector \(\boldsymbol{u}_{\mathrm{LOD}}\) associated with the corrected multiscale basis \(\{\Psi_H^{i}\}\). It first extracts fine-grid features by a convolutional encoder built from Conv–GroupNorm–GELU and residual blocks and then performs nonlinear mixing on the coarse LOD grid by residual layers, so that the output is restricted to the interior coarse degrees of freedom. Training is supervised by precomputed LOD coefficient vectors and uses the energy-based objective in 12 , which weights coefficient errors by the LOD stiffness matrix and therefore reflects their contribution to the reconstructed solution. At inference time, the predicted coefficients are combined with the exact corrected basis functions to reconstruct the solution, which is then evaluated on a uniform \(128\times128\) grid.

LOD-MSNO (Joint) combines the coefficient model with a learned basis-correction model, so that both the expansion coefficients \(\boldsymbol{u}_{\mathrm{LOD}}\) and the coefficient-dependent localized LOD basis \({\lbrace \Psi_H^j \rbrace}_{j=1}^{N_H}\) are predicted from the permeability field. The basis model is realized by PatchCorrectionNet, which acts locally on oversampling patches and predicts the correction coefficients associated with the three local coarse basis functions of each triangular coarse element. Because boundary and interior patches have different active fine degrees of freedom, each patch is embedded into a common rectangular patch canvas using padding, and the network receives the normalized log-permeability together with input and target support masks. Architecturally, PatchCorrectionNet is a U-Net-style convolutional encoder–decoder with skip connections, using the same Conv–GroupNorm–GELU and residual blocks as the coefficient model. The local patch predictions are assembled into a global approximation of the LOD correction operator and combined with the predicted coefficient vector to reconstruct the fine-scale solution.

All baseline models are trained directly on pairs \(\{(\kappa_i, u_h^{(i)})\}\) of permeability fields and corresponding fine-scale reference solutions computed by (3 ) on the fine grid with mesh size \(h\), using the relative \(L^2\) error as the training objective. For the models LOD-MSNO, FNO, CNO, and UNO, we apply the same input preprocessing and normalization pipeline, consisting of a logarithmic transformation of the permeability field, gradient-magnitude and coordinate channels, and channel-wise input as well as per-component output normalization (see Appendix 10.1 for details). For PI-DeepONet and Transolver, which do not support additional input channels, the gradient-magnitude and coordinate augmentations (add_grad and add_coords) are therefore omitted. For evaluation, we report the mean and standard deviation of the relative \(L^2\) error (Rel-\(L^2\)) computed on the fine grid, where \(u\) denotes the reference solution and \(\widehat{u}\) the model prediction.

\begin{table}[t]

Model	Lognormal1	Lognormal2	Quantile	Checkerboard
FNO [6]	\(0.049 \pm 0.001\)	\(0.163 \pm 0.002\)	\(0.254 \pm 0.010\)	\(0.064 \pm 0.000\)
CNO [19]	\(\underline{0.018 \pm 0.001}\)	\(0.093 \pm 0.003\)	\(\underline{0.120 \pm 0.007}\)	\(\mathbf{0.010 \pm 0.001}\)
UNO [20]	\(0.044 \pm 0.000\)	\(0.160 \pm 0.003\)	\(0.239 \pm 0.008\)	\(0.064 \pm 0.001\)
PI-DeepONet [17]	\(0.138 \pm 0.001\)	\(2.028 \pm 0.084\)	\(1.282 \pm 0.015\)	\(0.062 \pm 0.000\)
Transolver [8]	\(0.082 \pm 0.022\)	\(0.341 \pm 0.016\)	\(0.610 \pm 0.050\)	\(0.053 \pm 0.003\)
PINO [21]	\(0.025 \pm 0.000\)	\(0.145 \pm 0.002\)	\(0.203 \pm 0.012\)	\(0.035 \pm 0.000\)
FEONet [26]	\(0.099 \pm 0.000\)	\(1.805 \pm 0.073\)	\(0.129 \pm 0.020\)	\(0.668 \pm 0.109\)
LOD-MSNO (Coeff)	\(\mathbf{0.012 \pm 0.002}\)	\(\mathbf{0.042 \pm 0.006}\)	\(\mathbf{0.118 \pm 0.017}\)	\(\underline{0.018 \pm 0.003}\)
LOD-MSNO (Joint)	\(0.029 \pm 0.001\)	\(\underline{0.068 \pm 0.002}\)	\(0.302 \pm 0.016\)	\(0.029 \pm 0.001\)

\vskip

-0.1in \end{table}

The aggregated results are summarized in Table ¿tbl:tab:main95comparison?. The results indicate that LOD-MSNO(Coeff) attains competitive and, in three out of four cases, improved relative \(L^2\)-error accuracy compared to the considered baselines, with particularly strong performance on the multiscale lognormal datasets, while remaining comparable to the best-performing methods on the other coefficient configurations. In comparison to FEONet and PI-DeepONet, the results highlight the advantage of employing problem-adapted multiscale basis functions. While FEONet relies on a fixed coarse \(P1\)-FEM basis and PI-DeepONet learns a global set of basis functions via the trunk network, LOD-MSNO leverages LOD basis functions that are specifically tailored to the underlying coefficient field. This leads to consistently improved accuracy, particularly in strongly heterogeneous regimes. Moreover, compared to physics-informed approaches such as PINO and PI-DeepONet, the results suggest that incorporating the physical structure through a strong numerical prior is more effective than enforcing the PDE via autodifferentiation-based loss terms. In particular, constructing the solution using the LOD framework allows the model to directly encode relevant multiscale features of the problem, which appears to be a more robust way of injecting physical information than relying solely on PDE residual minimization. The results for LOD-MSNO (Joint) further show that the fully learned surrogate remains competitive despite also replacing the localized corrector computations. In particular, on the high-contrast lognormal2 dataset, LOD-MSNO (Joint) achieves a mean relative \(L^2\) error of \(0.068\), outperforming all benchmark models and demonstrating the potential of the full surrogate in high-contrast regimes.

5.2.2 Ablation study↩︎

\begin{table}[t]

3pt

Dataset	Training setting	Variants
3-6		LOD-MSNO (Coeff)	w/o log +normalization	w/o add_grad +add_coords	w/o all
	log + normalization
	add_grad + add_coords
Lognormal1		0.025	0.080	0.034	0.236
Lognormal2		0.042	0.135	0.096	1.215
Quantile		0.089	0.247	0.319	0.313
Checkerboard		0.020	0.020	0.056	0.188

\vskip

-0.1in \end{table}

We conduct an ablation study to assess the impact of the training loss function on the performance and generalization ability of LOD-MSNO (Coeff). In particular, we train the model using three alternative losses: the standard \(L^2\) loss, the relative \(L^2\) loss, and the \(L^1\) loss. We specifically use LOD-MSNO (Coeff) to isolate the effect of the choice of loss function from LOD-MSNO (Joint), in which the error predicting the basis also influences the relative L2 error on the fine grid. All models are trained on 500 samples with lognormal2 and quantile-type permeability fields, and we report the evolution of the test relative \(L^2\) error during training. As shown in Figure 5, training with the energy-based loss yields the lowest final error for both coefficient classes, while the \(L^1\) and (relative) \(L^2\) saturate at higher error levels. Especially in the lognormal2 case, the energy-loss exhibits the most stable convergence behavior, whereas the \(L^2\)- and \(L^1\)-based training runs display larger fluctuations and higher asymptotic errors. Table ¿tbl:tab:ablation95errors? reports the impact of the proposed preprocessing and channel augmentation strategies on the performance of LOD-MSNO (Coeff). Here, all models are trained using the energy-norm loss on 500 samples from the Quantile dataset. Across all datasets, the full configuration—combining logarithmic scaling, normalization, and the additional gradient-magnitude and coordinate channels—consistently yields the lowest relative \(L^2\) errors. Removing either the log-scaling and normalization or the auxiliary channels leads to a noticeable degradation in accuracy, while omitting all preprocessing steps results in a substantial loss of performance. These results demonstrate that both the normalization strategy and the inclusion of gradient and coordinate information play an important role in enabling accurate and stable learning for LOD-MSNO(Coeff).

5.3 Comparison of training and inference time↩︎

Table ¿tbl:tab:time95comparison? summarizes the data-generation costs for different mesh resolutions; all times are measured in seconds per sample. For the LOD basis construction, we report the local patch-corrector cost, the global basis assembly cost, and their sum parallelized over 16 CPU cores, the sequential global basis assembly cost, and their sum. The localized corrector problems \(Q_{T,h}^{\ell}\) are embarrassingly parallel—each depends only on the fine-grid stiffness matrix restricted to its patch—so the patch-corrector step distributes exactly across cores; scaling to more than 16 cores would reduce this component proportionally. By contrast, assembling the LOD basis (given by the rows of the matrix \(G = P_1 - Q_h^{\ell}\), which is the difference of the prolongation matrix P1 that encodes the coarse nodal basis functions and the correction operator \(Q_h^{\ell}\)) via global sparse-matrix stacking and multiplication across all patches such that is an aggregation step that does not decompose into independent patch-level tasks and is therefore measured without inter-task parallelization. Classical LOD timings were obtained on a 16-core CPU; all neural-network timings (training and inference) were measured on an NVIDIA A100 80GB GPU.

\begin{table}[t]

Resolution	Fine FEM	LOD Basis			LOD Coeff
3-5(lr)6-8		Local patch correctors	Global basis assembly	Total	Coeff. assembly	Coeff. solve	Total
\(h=7,\;H=4\)	0.1135	0.7842	0.4890	1.2732	0.2145	0.0020	0.2165
\(h=8,\;H=4\)	0.5590	4.2527	1.8327	6.0854	0.8500	0.0021	0.8521
\(h=8,\;H=5\)	0.5606	4.5842	4.7611	9.3453	1.2688	0.0752	1.3440
\(h=9,\;H=6\)	4.3828	20.0033	72.9008	92.9041	6.9178	3.1047	10.0225

\vskip

-0.1in \end{table}

The results show that the localized basis-corrector computations dominate the classical LOD data-generation cost, especially as the fine mesh is refined. With 16-core parallelization of the local corrector step, the per-sample patch-corrector time is substantially reduced compared to the sequential cost; using more than 16 cores would decrease this cost further proportionally. The global basis assembly and the coefficient steps remain sequential, and for fine coarse spaces, the global assembly itself becomes a dominant cost—at \(h=9\), \(H=6\), it accounts for the majority of the basis construction time. In particular, as \(H\) is refined, the coarse LOD solve becomes comparable in cost to the fine FEM solve despite involving fewer degrees of freedom, which highlights the effect of the non-sparsity of the LOD system as a major computational bottleneck. This observation motivates the coefficient surrogate, whose direct inference of the expansion coefficients avoids repeatedly assembling and solving this dense coarse system and thereby contributes directly to the speedup of LOD-MSNO. These timings motivate the two surrogate components of LOD-MSNO: the basis model targets the expensive localized corrector computations, whereas the coefficient model targets the assembly and solution of the coarse LOD system. Consequently, the timing results indicate that LOD-MSNO indeed addresses the main computational bottlenecks of the classical LOD method discussed in Section 3.3.

Table ¿tbl:tab:training95inference95comparison? reports the corresponding neural-network training and single-sample inference times for LOD-MSNO and CNO. All neural-network models were trained using 1000 training samples, and the coefficient model and CNO were trained for the same number of epochs. We compare against CNO because it is the strongest benchmark model in terms of accuracy among the standard neural-operator baselines considered in Table ¿tbl:tab:main95comparison?.

\begin{table}[t]

3pt

Resolution	Training time			Inference
2-4(lr)5-7	Coeff	Basis	CNO	Coeff	Basis	CNO
\(h=7,\;H=4\)	1186	4409	10366	0.0036	0.0347	0.0061
\(h=8,\;H=4\)	3922	12290	36076	0.0036	0.9686	0.0064
\(h=8,\;H=5\)	4508	17652	36040	0.0036	1.3884	0.0064
\(h=9,\;H=6\)	27334	71862	137714	0.0086	5.5559	0.0207

\vskip

-0.1in \end{table}

The results show that both LOD-MSNO components combined are faster to train than CNO across all tested resolutions. At inference time, the neural surrogates are much faster than the corresponding traditional LOD computations, as can be inferred from the per-sample LOD basis and coefficient generation times in Table ¿tbl:tab:time95comparison?. This speedup stems from replacing the repeated localized corrector solves and coarse LOD solves with neural-network predictions, and holds even when the classical pipeline is parallelized over 16 CPU cores for the corrector step. Consequently, the inference-time results indicate that LOD-MSNO is capable of addressing the main computational bottlenecks of the classical LOD method discussed in Section 3.3. The CNO inference time is lower than that of the basis model because CNO directly outputs the solution as an image on the grid. By contrast, the basis model predicts local basis corrections patch by patch and must assemble the global correction operator and reconstructed solution from these individual local predictions.

6 Conclusion and outlook↩︎

In this work, we introduced LOD-MSNO, a hybrid neural-operator framework for multiscale Darcy-Flow Problems that combines the problem-adapted approximation structure of the Localized Orthogonal Decomposition method with deep learning-based surrogate modeling. The method is designed based on the Localized Orthogonal Decomposition (LOD) method and aims to replace the expensive numerical computations associated with it through a mainly data-driven approach. In the coefficient-only variant, LOD-MSNO (Coeff) predicts the LOD expansion coefficients and reconstructs the fine-scale solution using the exact corrected multiscale basis, thereby avoiding the assembly and solution of a non-sparse system of equations. In the joint variant, LOD-MSNO (Joint) additionally learns the localized basis corrections, thereby also targeting the second computational bottleneck of the LOD method: the coarse LOD solve and the construction of localized correctors. Through a rigorous error analysis, we gave a precise convergence bound for our models. Numerical experiments also show that using the LOD method as a modeling prior is beneficial for challenging rough and high-contrast permeability fields for the Darcy-flow equation. The comparison with other physics-informed baselines further indicates that explicitly incorporating problem-adapted LOD basis information can be more effective than relying only on PDE residuals for physically consistent solutions.

The following directions remain open for future work. A stronger theoretical analysis of the basis-learning component would be valuable, especially for quantifying how errors in the learned local correction operators propagate to the assembled multiscale basis and the final fine-scale reconstruction. Another important direction is to improve the efficiency of the joint model during inference, for example by reducing patch-wise assembly overhead or by exploiting parallel structure in the local corrector predictions. Lastly, since LOD methods are available for a wide range of problem classes, including elliptic eigenvalue problems and parabolic equations, extending LOD-MSNO to these settings constitutes another promising avenue for future work.

Acknowledgements↩︎

This work was supported by Samsung Electronics Co., Ltd (IO250307-12279-01), by the Institute of Information & Communications Technology Planning & Evaluation (IITP) grant funded by the Korea government (MSIT) [RS-2021-II211341, Artificial Intelligence Graduate School Program (Chung-Ang University)], and by the Chung-Ang University Young Scientist Scholarship in 2025.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process↩︎

During the preparation of this work, the author(s) used Prism AI in order to polish and improve the scientific language and readability of human-written texts, as well as to assist with the formatting of formulas and tables. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the published article.

7 Theoretical Analysis↩︎

7.1 SPD-property of the LOD-stiffness matrix↩︎

7.2 Proof of Theorem 2 - Approximation error for the energy-based loss function↩︎

To prove Theorem 2, we shall use the following Lemma:

7.3 Proof of Theorem 3 - Generalization error for the energy-based loss function↩︎

8 Baseline models↩︎

8.1 Common input–output format and training losses↩︎

All baselines learn the solution operator for (1 ) on a regular \(s\times s\) grid. The input is a coefficient field (and optional engineered channels), and the output is the scalar solution field. Concretely, the network input is a tensor in \(\mathbb{R}^{B\times s\times s\times C_{\mathrm{in}}}\) (or equivalently \(B\times C_{\mathrm{in}}\times s\times s\)), and the output prediction is \(u_\theta \in \mathbb{R}^{B\times s\times s}\). Here, \(C_{\mathrm{in}}\) is the number of input channels used in the experiment. All baselines except FEONet are trained on fine-mesh finite element reference solutions \(u_h\) defined by 3 , i.e \(s=\frac{1}{h}\). FEONet is treated separately because it is trained through its finite-element weak-form residual rather than by direct supervision on these fine-mesh solution labels.

For supervised training, we use the standard data loss \[\mathcal{L}_{\mathrm{data}}(\theta)=\|u_\theta-u\|_2^2.\] For physics-informed training of PI-DeepONet and PINO, we additionally enforce the Darcy equation \[-\nabla \cdot \big(\kappa(x)\nabla u(x)\big)=f(x)\] through coordinate-based automatic differentiation rather than a finite-volume/finite-difference stencil. Let \[\mathcal{G}=\{x_{ij}\}_{i,j=1}^{s}\] denote the regular grid points, and let \(u_\theta(x;\kappa)\) be the network prediction evaluated at coordinate \(x\in\Omega\). We define the pointwise residual by \[r_\theta(x;\kappa) = -\partial_{x_1}\!\left(\kappa(x)\,\partial_{x_1}u_\theta(x;\kappa)\right) -\partial_{x_2}\!\left(\kappa(x)\,\partial_{x_2}u_\theta(x;\kappa)\right) -f(x).\] In the implementation, the spatial coordinates are treated as differentiable inputs, \(\nabla u_\theta\) is obtained by automatic differentiation with respect to those coordinates, and \(\kappa(x)\) is evaluated from the raw coefficient field at the same coordinates by differentiable interpolation. The physics-informed loss functions are then \[\mathcal{L}_{\mathrm{pde}}(\theta) = \frac{1}{|\mathcal{G}|}\sum_{x\in\mathcal{G}} |r_\theta(x;\kappa)|^2, \qquad \mathcal{L}_{\mathrm{PI}}(\theta) = \mathcal{L}_{\mathrm{data}}(\theta)+\lambda \mathcal{L}_{\mathrm{pde}}(\theta).\]

8.2 Fourier neural operator (FNO)↩︎

The Fourier Neural Operator (FNO), introduced by Li et al. [6], is a neural operator architecture for learning mappings between infinite-dimensional function spaces, with a particular focus on approximating solution operators of parametric partial differential equations. The key idea of the FNO is to parameterize an integral operator in Fourier space, enabling efficient global convolution and resolution-invariant generalization.

Given a function \(v_t\) defined on a spatial domain, the FNO constructs a learnable integral operator \(\mathcal{K}(\phi)\) of the form \[\left( \mathcal{K}(\phi)v_t \right)(x) = \mathcal{F}^{-1} \left( \mathcal{F}(\kappa_{\phi}) \cdot (\mathcal{F}v_t) \right)(x),\] where \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) denote the Fourier transform and its inverse, respectively, and \(\kappa_{\phi}\) is a learnable kernel. This formulation corresponds to a global convolution in physical space and allows the operator to capture long-range dependencies efficiently.

In practice, the continuous Fourier transforms are replaced by discrete fast Fourier transforms. Specifically:

• The operators \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) are replaced by the discrete fast Fourier transforms \(\hat{\mathcal{F}}\) and \(\hat{\mathcal{F}}^{-1}\).

• The input \(\hat{v} \in \mathbb{R}^{s_1 \times \dots \times s_d \times n}\) represents the evaluation of the function \(v\) on a uniform grid with \(s_j \in \mathbb{N}\) points in each spatial dimension.

• The Fourier transform of the kernel \(\mathcal{F}(\kappa)\) is replaced by a complex-valued weight tensor \(\boldsymbol{T} \in \mathbb{C}^{s_1 \times \dots \times s_d \times n \times n}\), whose entries correspond to the learnable Fourier modes of the kernel.

To reduce computational complexity and improve generalization, the FNO typically employs a low-frequency truncation in Fourier space. Concretely, only a prescribed number of low-frequency modes are retained in each spatial dimension, while the remaining high-frequency modes are set to zero. With this truncation, one Fourier layer of the FNO is given by \[\hat{v} \;\mapsto\; \hat{\mathcal{F}}^{-1} \left( \boldsymbol{T}_{\text{low}} \cdot \hat{\mathcal{F}}(\hat{v}) \right),\] where \(\boldsymbol{T}_{\text{low}}\) denotes the restriction of \(\boldsymbol{T}\) to the retained low-frequency Fourier modes. By stacking multiple such Fourier layers and interleaving them with pointwise nonlinearities, the Fourier Neural Operator is able to approximate nonlinear operators with global receptive fields while remaining computationally efficient.

8.2.0.1 Experimental configuration

In our experiments, the FNO uses \(\textcolor{black}{4}\) Fourier layers with channel width \(\textcolor{black}{64}\), retaining \(\textcolor{black}{12}\) low-frequency modes per spatial dimension in each spectral convolution layer.

8.3 Convolutional neural operator (CNO)↩︎

The Convolutional Neural Operator (CNO) [19] is designed to approximate operators \(\mathbf{G} : \mathcal{B}_w(D,\mathbb{R}^{d_X}) \to \mathcal{B}_w(D,\mathbb{R}^{d_Y})\) between band-limited function spaces. By restricting both the input and output to band-limited spaces, the CNO mitigates aliasing errors and enables stable learning across different spatial resolutions.

The CNO represents the target operator as a composition of lifting, convolutional operator layers, and projection, expressed as \[\mathcal{G} : u \mapsto P(u) = v_0 \mapsto v_1 \mapsto \cdots \mapsto v_L \mapsto Q(v_L) = \bar{u},\] where \(P\) is a lifting operator that maps the input function \(u\) to a higher-dimensional feature representation, \(Q\) is a projection operator to the target space, and the intermediate states \(\{v_\ell\}_{\ell=0}^L\) are defined recursively by \[v_{\ell+1} = \mathcal{P}_{\ell} \circ \Sigma_{\ell} \circ \mathcal{K}_{\ell}(v_{\ell}), \qquad 1 \le \ell \le L-1.\] Here, \(\mathcal{K}_{\ell}\) denotes a convolutional operator, \(\Sigma_{\ell}\) a pointwise nonlinearity, and \(\mathcal{P}_{\ell}\) a resolution-changing operator.

Architecturally, the CNO adopts a modified operator U-Net structure composed of an encoder that expands the feature width and contracts the spatial height, enabling multiscale feature extraction, and a decoder that reconstructs high-resolution outputs by successively upsampling the feature maps while reducing the channel dimension. To compensate for the potential loss of high-frequency information caused by upsampling, feature maps from earlier encoder stages are concatenated with decoder features, ensuring the preservation of fine-scale details.

8.3.0.1 Experimental configuration

We use a \(\textcolor{black}{4}\)-level encoder–decoder CNO with base channel width \(\textcolor{black}{64}\) and \(\textcolor{black}{2}\) residual (R-)blocks at each resolution.

8.4 U-shaped neural operator (UNO)↩︎

UNO is a Fourier-based neural operator architecture that extends the Fourier Neural Operator (FNO) with a U-Net-style encoder-decoder structure and skip connections. UNO progressively contracts and expands the spatial domain of intermediate functions, enabling deeper and more memory-efficient operator learning while preserving global interactions through spectral integral operators. Its core operations are performed in the Fourier domain, making it particularly effective for learning global solution operators of PDEs.

In contrast, CNO replaces Fourier integral operators with purely convolutional operators defined in physical space. CNOs rely on hierarchical convolution, interpolation, and local kernel operations to achieve resolution-invariant operator learning. Unlike UNO, CNO does not explicitly use spectral transforms or global Fourier modes; instead, it emphasizes locality, translation equivariance, and scalability through convolutional architectures similar to CNNs.

Therefore, the key distinction is that UNO learns operators primarily through global spectral interactions using Fourier operators and domain contraction/expansion, whereas CNO learns operators through local multiscale convolutions directly in physical space.

8.4.0.1 Experimental configuration

We use an \(\textcolor{black}{4}\)-stage UNO with base channel width \(\textcolor{black}{64}\) in the lifting layer and scale the width across encoder/decoder stages accordingly.

8.5 Deep operator network (DeepONet / PI-DeepONet)↩︎

Deep Operator Networks (DeepONets) [35] learn nonlinear operators by combining two subnetworks: a branch net that encodes the input function and a trunk net that encodes the query location. For an input coefficient field \(\kappa\) and a spatial point \(y\in\Omega\), DeepONet represents the solution operator in the form \[u_\theta(y;\kappa) = \sum_{j=1}^{p} b_{\theta,j}(\kappa)\, t_{\theta,j}(y), \label{eq:deeponet95separable}\tag{26}\] where \(b_{\theta}(\kappa)\in\mathbb{R}^{p}\) is produced by the branch net and \(t_{\theta}(y)\in\mathbb{R}^{p}\) is produced by the trunk net. Evaluating 26 over all grid points yields the full predicted field.

Physics-informed variants of DeepONet (often referred to as PI-DeepONet or PI-DON) incorporate the governing PDE as an additional training signal. In our experiments, we adopt a PI-DeepONet-style training of DeepONet by using the common PI objective in Section 8.1, where the physics loss is computed by automatic differentiation with respect to the trunk coordinates.

8.5.0.1 Experimental configuration

For the plain DeepONet baseline, we set the latent dimension to \(p=\textcolor{black}{128}\). Both the branch and trunk nets are MLPs of depth \(\textcolor{black}{4}\) and \(\textcolor{black}{4}\), respectively, with hidden width \(\textcolor{black}{128}\).

For the PI-DeepONet baseline, we use branch layers \((d_b,100,100,100,100)\) and trunk layers \((2,100,100,100,100)\), where \(d_b=C_{\mathrm{in}}s^2\) is the flattened branch input dimension.

8.6 Transolver↩︎

Transolver [8] is an operator learning architecture built on physics-attention, which compresses dense attention over grid tokens into a two-stage slicing–token mechanism. Let \(x\in\mathbb{R}^{N\times d}\) denote the per-point features on a grid (\(N=s^2\)). Transolver first assigns each grid point to one of \(S\) latent slices by \[A=\mathrm{softmax}\!\left(\frac{W_s x}{\tau}\right)\in\mathbb{R}^{N\times S},\] then forms slice tokens by normalized pooling, \[T=\frac{A^\top x}{A^\top \mathbf{1}}\in\mathbb{R}^{S\times d}.\] Multi-head self-attention is applied among the \(S\) tokens (cost \(\mathcal{O}(S^2)\)), and the updated tokens are projected back to the grid by deslicing, \(\tilde{x}=A T'\). This reduces the quadratic cost in \(N\) to \(\mathcal{O}(NS+S^2)\) while retaining global information flow, which is beneficial for elliptic problems where long-range dependencies are intrinsic.

8.6.0.1 Experimental configuration

We use feature dimension \(d=\textcolor{black}{128}\), number of slices \(S=\textcolor{black}{64}\), number of attention heads \(\textcolor{black}{8}\), and temperature \(\tau=\textcolor{black}{0.5}\) for the slicing softmax.

8.7 Physics-informed Neural Operator (PINO)↩︎

The Physics-Informed Neural Operator (PINO) combines an operator-learning backbone with a PDE residual during training. In our implementation, the PINO baseline uses an FNO-style 2D architecture on the regular \(s\times s\) grid. Concretely, if \(v_\ell\) denotes the hidden feature map at layer \(\ell\), one PINO block takes the form \[v_{\ell+1}(x) = \sigma\!\left(\mathcal{K}_\ell(v_\ell)(x)+W_\ell v_\ell(x)\right),\] where \(\mathcal{K}_\ell\) is a Fourier spectral convolution retaining only low-frequency modes and \(W_\ell\) is a learnable pointwise linear map. Stacking such blocks yields a neural operator with a global receptive field and efficient spectral mixing.

For physics-informed training, PINO uses the same PI objective from Section 8.1.

8.7.0.1 Experimental configuration

We use four spectral blocks, hidden width \(64\), retain \(12\) Fourier modes in each spatial dimension at every block, use a pointwise projection dimension of \(128\), GELU activations, and no domain padding.

8.8 Finite-Element Operator Network (FEONet)↩︎

FEONet is a neural operator model grounded in the finite element method (FEM), using nodal basis functions \(\{\phi_i\}\) to approximate solutions as a linear combination \[u = \sum_{i=1}^{N_H} \alpha_i \phi_i\] where the index \(i\) denotes the nodes of a triangulation of the physical domain and \(N_H\) the number of (coarse) nodes. The basic idea of FEONet is to predict the FEM-coefficients \(\{\alpha_i\}\). For the Darcy equation 1 , the training procedure of FEONet is based on the weak formulation, which reads: Find \(u \in V\) such that \[\int_\Omega \kappa \nabla u \cdot \nabla v \, dx = \int_\Omega f v \, dx \quad \forall v \in V.\] Using the FEM ansatz \(u = \sum_j \alpha_j \phi_j\) and testing with \(\phi_i\), one obtains the linear system \(A \alpha = b\) with the stiffness matrix \(A = (A_{ij})\) and the load vector \(f = (f_i)\) defined as \[A_{ij} = \int_\Omega \kappa \nabla \phi_j \cdot \nabla \phi_i \, dx, \quad b_i = \int_\Omega f \phi_i \, dx.\] for \(i,j=1,....,N_H\). FEONet is then trained by sampling permeability fields \(\kappa_i\) and computing corresponding stiffness and load vector samples \(\{(A_i(\kappa_i), b_i(\kappa_i))\}_{i=1}^N\) and then using the weak-form residual loss over samples with index \(i\): \[\mathcal{L} = \frac{1}{N} \sum_{i=1}^N \| A_i \alpha_i - b_i \|.\] Note that, by using this loss function, we do not need the ground-truth FEM-coefficients \(\alpha_i\) to train our model. This way, FEONet can be interpreted as a hybrid Operator Learning model that combines neural networks with the conventional FEM-method, which is data-free and fully physics-informed.

8.8.0.1 Experimental configuration

In our implementation of FEONet, we use a CNN-based feature extractor for the input \(\kappa\), followed by a Multi-Layer Perceptron to predict the coefficients \(\alpha_i\) based on these features. For comparability with LOD-MSNO, the input of the model is the permeability field \(\kappa\) on a grid with the fine resolution \(h\) and the output is the coefficients \(\alpha_i\) for the coarse-scale-basis functions \(\phi_i\), \(i=1,...,N_H\).

9 Datasets↩︎

9.0.0.1 Coefficient families

We provide more detailed information about the used datasets here:

• Quantile and filtered quantile coefficients. We first sample a smooth Gaussian random field \(g\) on the fine nodes using the squared-exponential covariance \[\mathbb{E}\bigl[g(x)g(y)\bigr] = \sigma_{\mathrm{SE}}^2 \exp\!\left(-\frac{\|x-y\|_2^2}{2\ell_{\mathrm{SE}}^2}\right),\] where \(\ell_{\mathrm{SE}}=\texttt{se\_length}\) and \(\sigma_{\mathrm{SE}}=\texttt{se\_sigma}\). The field is exponentiated, \(a=\exp(g)\), and empirical quantile bins of \(a\) are mapped to a randomly shuffled set of prescribed permeability values \[\{0.01,4,8,12,66,100\}\] by default. The filtered_quantile variant additionally blends the nodal values toward the background value \(4\) near the boundary using a smooth boundary weight. This is the default coefficient type in the generation script.

• Lognormal Matérn coefficients. For coeff_type=lognormal, we sample a Matérn Gaussian field \(Z\) on the fine nodes and define \(\kappa=\exp(Z)\). The covariance is \[\mathbb{E}\bigl[Z(x)Z(x')\bigr] = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)} \left( \frac{\sqrt{2\nu}}{\ell}\|x-x'\|_2 \right)^{\nu} K_{\nu}\!\left( \frac{\sqrt{2\nu}}{\ell}\|x-x'\|_2 \right),\] where \(K_{\nu}\) denotes the modified Bessel function of the second kind. The default Matérn parameters are \(\nu=0.5\), correlation length \(\ell=0.3\), variance parameter \(\sigma=2.0\), and jitter \(10^{-8}\), with the same implementation also allowing the lower-contrast lognormal setting by changing these arguments.

• Fine checkerboard coefficients. For coeff_type=fine_checkerboard, values are sampled independently from the prescribed permeability levels on the fine \(h\)-scale grid cells. Nodal coefficient values are obtained by indexing into the corresponding fine cell.

9.0.0.2 Per-sample LOD assembly

For each generated coefficient realization, the script performs the following steps:

1. Fine stiffness assembly. The elementwise coefficient is used to assemble a discontinuous Galerkin element stiffness matrix, which is then mapped to the continuous fine space to obtain the fine stiffness matrix \(\boldsymbol{A}_h\). The load corresponds to the constant source \(f\equiv 1\) and is assembled using the fine mass matrix.

2. Localized correctors. For each coarse element \(T\), the localized corrector is computed on the precomputed fine patch degrees of freedom. In matrix form, the implementation solves the constrained local corrector problem using the local stiffness block, the local right-hand side induced by the three coarse basis functions on \(T\), and the patch interpolation constraint \(\boldsymbol{I}_{H,T}\). The resulting local corrector matrix has shape \(N_{\ell,h}(T)\times 3\) and is stored in the flattened array corr_values, with corr_ptr marking the offsets of each patch.

3. Global multiscale operator. The local correctors are inserted into global fine-grid columns and assembled into the localized correction matrix \(\boldsymbol{C}^{\ell}\). The multiscale reconstruction matrix is then \[\boldsymbol{G}=\boldsymbol{P}_1-\boldsymbol{C}^{\ell}.\] After restricting to interior coarse degrees of freedom, the LOD stiffness matrix and load vector are formed as \[\boldsymbol{A}_{\mathrm{LOD}} = \boldsymbol{G}_0^{\top}\boldsymbol{A}_h\boldsymbol{G}_0, \qquad \boldsymbol{f}_{\mathrm{LOD}} = \boldsymbol{G}_0^{\top}\boldsymbol{b}_h.\]

4. Coarse solve and fine reconstructions. The interior LOD system is solved for \(\boldsymbol{u}_{\mathrm{LOD}}\), the full coarse vector is padded with zero boundary values, and the fine-grid multiscale solution is reconstructed as \[\boldsymbol{u}_{\mathrm{ms}}=\boldsymbol{G}\boldsymbol{u}_{H}.\] For diagnostic and evaluation purposes, the fine FEM reference solution \(\boldsymbol{u}_h\) is also computed by solving the fine interior system with homogeneous Dirichlet boundary conditions.

10 Implementation details↩︎

10.1 Neural architecture for LOD-MSNO (Coeff)↩︎

10.1.0.1 Input preprocessing

We form a multi-channel tensor \(X\in\mathbb{R}^{B\times C_{\mathrm{in}}\times n_h\times n_h}\) from the raw field \(a_h\). Denote by \(b\) the base channel: \[b(i, j) = \log(\max\{a_h(i,j),\varepsilon\})\] where \(\varepsilon>0\) is a small constant to avoid numerical issues. The log transform compresses the dynamic range of high-contrast permeability fields. This is consistent with the common modeling practice in Darcy-flow settings [36].

10.1.0.2 Gradient-magnitude channel (add_grad).

To expose sharp interfaces/discontinuities in \(a_h\), we append a gradient-magnitude channel computed by forward differences (with zero padding at the boundary): \[\begin{align} g_x(i,j) &= \begin{cases} b(i,j+1)-b(i,j), & j<n_h-1,\\ 0, & j=n_h-1, \end{cases} \\ g_y(i,j) &= \begin{cases} b(i+1,j)-b(i,j), & i<n_h-1,\\ 0, & i=n_h-1, \end{cases} \\ g(i,j) &= \sqrt{g_x(i,j)^2 + g_y(i,j)^2 + \delta}, \label{eq:preproc95gradmag} \end{align}\tag{27}\] with \(\delta>0\) small.

10.1.0.3 Coordinate channels (add_coords).

Convolutional networks are translation-equivariant, while boundary-value problems on bounded domains are position-dependent. We therefore append explicit coordinates (CoordConv-style positional encoding [37]): \[x(j)=-1+\frac{2j}{n_h-1}, \qquad y(i)=-1+\frac{2i}{n_h-1}, \label{eq:preproc95coords}\tag{28}\]

10.1.0.4 Channel count

The number of input channels is \[C_{\mathrm{in}} \;=\; 1 \;+\; \mathbb{I}[\texttt{add\_grad}] \;+\; 2\,\mathbb{I}[\texttt{add\_coords}], \label{eq:Cin}\tag{29}\] and the final preprocessed tensor is formed by concatenation along the channel dimension.

Let \(X\in\mathbb{R}^{B\times C_{\mathrm{in}}\times n_h\times n_h}\) be the preprocessed input.

10.1.0.5 Channel-wise input normalization.

We compute per-channel statistics \((\mu_c,\sigma_c)\) on the training set after preprocessing: \[\mu_c \;=\; \mathbb{E}[X_{c}],\qquad \sigma_c^2 \;=\; \mathbb{E}[X_{c}^2]-\mu_c^2, \label{eq:x95stats}\tag{30}\] where expectations are taken over all training samples and spatial locations. We then normalize as \[\widetilde{X}_{b,c,i,j} \;=\; \frac{X_{b,c,i,j}-\mu_c}{\sigma_c+\epsilon}. \label{eq:x95norm}\tag{31}\] The values \(\mu,\sigma\) are stored as non-trainable buffers and reused at test time.

10.1.0.6 Per-component output normalization.

Let \(u_{\mathrm{int}}\in\mathbb{R}^{D}\) denote the target interior coefficient vector. We compute per-component statistics \((\mu^{(u)}_k,\sigma^{(u)}_k)\) on the training set and normalize: \[\widetilde{u}_k \;=\; \frac{u_{\mathrm{int},k}-\mu^{(u)}_k}{\sigma^{(u)}_k+\epsilon}. \label{eq:u95norm}\tag{32}\] The core network is trained to predict \(\widehat{\widetilde{u}}\); inference returns the decoded coefficients \(\widehat{u}_{\mathrm{int}}\) via the inverse transform of 32 .

10.1.0.7 Building block for CNN-based model: GNAct.

We use the block \[\mathrm{GNAct}(z) \;=\; \mathrm{GELU}\!\left(\mathrm{GN}\!\left(\mathrm{Conv}_{k\times k,s}(z)\right)\right), \label{eq:gnact}\tag{33}\] where \(\mathrm{GN}\) denotes Group Normalization [38]. For GroupNorm, we use \(G=\min(8,C)\) groups for a \(C\)-channel tensor.

10.1.0.8 Residual block (ResBlock2D).

At a fixed resolution and channel width \(C\), we use a ResNet-style block: \[\mathrm{ResBlock}(x) = \mathrm{GELU}\left(x + \mathrm{GN} \left( \mathrm{Conv}_{3 \times 3, 1}\left( \mathrm{GNAct}(x) \right) \right) \right) \label{eq:resblock}\tag{34}\]

10.1.0.9 Hyperparameters used.

We use \(\texttt{base}=32\), GroupNorm with \(G=\min(8,C)\) groups, \(\texttt{add\_grad}=\texttt{True}\), and \(\texttt{add\_coords}=\texttt{True}\).

10.2 Neural architecture for predicting the basis correction↩︎

10.2.0.1 Patch-wise supervised learning target

Each training example corresponds to one coefficient realization and one coarse patch. For a fixed patch \(T\in\mathcal{T}_H\), the input is the coefficient field restricted to the localized region \(\mathcal{N}^{\ell}(T)\) and stored on a rectangular patch grid. Since boundary and interior patches can have different physical sizes and different active fine degrees of freedom, each patch is embedded into a common rectangular patch canvas using padding. Thus, all patches have the same input shape and can be processed by one convolutional network.

The target consists of the local correction coefficients \[\boldsymbol{w}_{T,j}(\omega) = \bigl(w_{T,j,1}(\omega),\ldots,w_{T,j,N_{\ell,h}(T)}(\omega)\bigr)^{\top}, \qquad j=1,\ldots,N_H,\] restricted in practice to the coarse basis functions whose support intersects \(T\). Since a two-dimensional triangular coarse element has three local coarse vertices, the implementation stores the local corrector target as an array with three output components per fine patch degree of freedom. Thus, for a patch with \(N_{\ell,h}(T)\) active fine degrees of freedom, the supervised target has shape \[Y_T(\omega)\in\mathbb{R}^{N_{\ell,h}(T)\times 3}.\] The basis dataset therefore consists of patch-sample pairs \[\mathcal{D}_{\mathrm{basis}} = \left\{ \left(\kappa_i|_{\mathcal{N}^{\ell}(T)},Y_T(\omega_i)\right) :\; i=1,\ldots,M,\; T\in\mathcal{T}_H \right\}.\] In the implementation, this dataset is flattened over samples and patches, so that one training epoch iterates over all pairs \((i,T)\).

10.2.0.2 PatchCorrectionNet architecture

The core of PatchCorrectionNet is a U-Net-style convolutional encoder–decoder with skip connections. It uses the same basic Conv–GroupNorm–GELU and residual blocks as the coefficient architecture, but preserves the patch-grid structure and returns a dense three-channel patch image. Starting from \[X_T\in\mathbb{R}^{B\times 3\times H_p\times W_p},\] the network applies the following stages:

• Stem: a non-downsampling DownStage mapping \(3\rightarrow c_1\) channels, followed by residual blocks.

• Encoder: three stride-two downsampling stages with channel widths \[\begin{align} c_1&=\texttt{base}, & c_2&=\min(2c_1,\texttt{max}), \\ c_3&=\min(4c_1,\texttt{max}), & c_4&=\min(8c_1,\texttt{max}). \end{align}\] Each stage consists of a Conv–GroupNorm–GELU block followed by blocks_per_stage residual blocks.

• Bottleneck: bottleneck_blocks residual blocks at width \(c_4\).

• Decoder: three upsampling stages using bilinear interpolation, concatenation with the corresponding encoder skip feature, a Conv–GroupNorm–GELU projection, and residual blocks.

• Head: a Conv–GroupNorm–GELU block, one residual block, and a final \(1\times1\) convolution producing \[\widehat{Y}_T^{\mathrm{grid}}(\omega) \in \mathbb{R}^{B\times 3\times H_p\times W_p}.\]

The three output channels correspond to the correction coefficients associated with the three local coarse basis functions on the patch. The architecture is fully convolutional on the patch canvas, and the use of skip connections is important because local correctors are spatially localized and may contain sharp features near coefficient jumps or patch boundaries.

10.2.0.3 Patch readout and constraint projection

The network output is a dense patch image, whereas the correction target is stored only on active fine degrees of freedom. Let \[\mathcal{C}_T=\{(x_m,y_m)\}_{m=1}^{N_{\ell,h}(T)}\] denote the list of local grid coordinates corresponding to these active fine degrees of freedom. The dense prediction is read out as \[\widehat{Y}_T(\omega) = \left( \widehat{Y}_T^{\mathrm{grid}}(:,y_m,x_m) \right)_{m=1}^{N_{\ell,h}(T)} \in\mathbb{R}^{N_{\ell,h}(T)\times 3}.\] Because the LOD correctors must lie in the local fine-scale space \(W^{\ell}(T)=\ker(\mathcal{I}_H)\), the implementation applies a hard algebraic projection before computing the loss. With the precomputed patch matrices denoted by \(\boldsymbol{P}_{T}^{\mathrm{left}}\) and \(\boldsymbol{P}_{T}^{\mathrm{right}}\), the projected prediction is \[\widehat{Y}_{T}^{\mathrm{proj}} = \widehat{Y}_T - \boldsymbol{P}_{T}^{\mathrm{left}} \boldsymbol{P}_{T}^{\mathrm{right}} \widehat{Y}_T,\] and the same projection is applied to the target \(Y_T\). This projection enforces compatibility with the kernel constraint induced by \(\mathcal{I}_H\) and prevents the network from being penalized for components outside the admissible local detail space.

10.2.0.4 Hyperparameters used

For PatchCorrectionNet, we use \[\begin{align} \texttt{base\_channels} &= 48, & \texttt{max\_channels} &= 256, & \texttt{blocks\_per\_stage} &= 2, \\ \texttt{bottleneck\_blocks} &= 4, & \texttt{dropout} &= 0.05. \end{align}\] The output dimension is fixed to three channels, corresponding to the three local coarse basis functions associated with a triangular coarse element. During validation, predictions from all patches are assembled into the global correction matrix \(\widehat{\boldsymbol{Q}}_h^{\ell}\) and combined with the coarse coefficients to reconstruct the fine-scale LOD solution.

References↩︎