A Generative Residual Enhanced Neural Operator Based on the Boundary Element Method for Accurate Metasurface Parameter Analysis

1  Introduction

A metasurface is an artificially designed metamaterial with subwavelength dimensions that follows the generalized Snell’s law and can achieve anomalous properties unattainable by traditional materials [1–4]. However, conventional material discovery and design processes are largely based on expensive, laborious, and time-consuming trial-and-error methods. The outcomes of these iterative processes not only rely heavily on human intuition and domain expertise, but also require a degree of luck [5,6]. Consequently, the ultimate goal of modern engineering is to rapidly design and analyze materials with outstanding performance [7,8]. However, whether it is material design or parameter analysis, the underlying physical principles involve complex multi-scale interactions, such as multi-scale phononic crystal design [9,10]. Therefore, parameter analysis of metasurfaces constitutes a challenging cross-scale problem [11,12].

In recent years, neural operator methods have attracted considerable attention due to their ability to learn mappings between infinite-dimensional function spaces, offering a promising alternative to traditional discretization-based solvers [13–16]. Among these methods, DeepONet [17,18] and the Fourier Neural Operator [19–21] (FNO) have demonstrated outstanding potential in approximating partial differential equation solvers in a wide range of scientific and engineering applications. However, existing neural operator models typically rely on large amounts of high-quality training data to achieve strong predictive performance [22,23]. For acoustic scattering problems, generating datasets requires extensive computation, especially when dealing with complex geometries and broadband metasurfaces [24]. These problems typically involve unbounded domains, which further increases the difficulty of dataset construction [25]. Therefore, achieving high-accuracy cross-scale prediction with reduced datasets remains a pressing challenge [26].

The Boundary Element Method (BEM) is an efficient numerical technique in which the core idea is to reduce the dimensionality of a problem by employing Boundary Integral Equations (BIEs)—for instance, transforming a three-dimensional (3D) problem into a two-dimensional (2D) boundary integral, or a two-dimensional problem into a one-dimensional boundary integral [27–29]. As a result, BEM requires discretization only on the boundary of the computational domain, which significantly reduces the computational cost. BEM is particularly well-suited for unbounded domain problems, such as computing the scattered acoustic fields of metasurfaces [30,31]. However, the fully populated dense system matrices generated by conventional BEM still impose high computational and storage costs. Moreover, large-scale analyses in varying parameter settings require repeated BEM computations, resulting in substantial computational resources [32].

In this work, we aim to leverage the dimensionality reduction feature of BEM to develop a high-accuracy neural operator for metasurface parameter analysis. In recent years, generative modeling methods—represented by Generative Adversarial Networks (GANs)—have demonstrated excellent performance in metasurface design [33–35]. Generative models are capable of explicitly capturing subtle details in high-dimensional data, which is often difficult for traditional discriminative models to achieve [36]. This inspired us to employ a generative mechanism to capture higher-order statistical features in the error space, thereby generating residual fields to compensate for the systematic bias in neural operator predictions. We refer to this framework as the Generative Residual Enhanced Neural Operator (GRE-NO). This strategy addresses the limitations of end-to-end fitting approaches in modeling complex error patterns under small boundary data scenarios, and it also offers a new path toward more accurate and robust neural operator predictions.

To the best of our knowledge, no existing work has explored such a generative residual-enhancement mechanism for boundary-based metasurface parameter analysis. Compared with existing methods that rely on large amounts of training data or fixed physical priors, GRE-NO offers a unified and efficient approach to metasurface parameter analysis with small-size data requirements and greater flexibility, while keeping the computational overhead manageable. Our method establishes an effective mapping from metasurface geometric parameters and incident wave parameters to boundary scattering field solutions, enabling accurate cross-scale analysis even under severe data scarcity. We validate the proposed approach through comprehensive numerical experiments on both two-dimensional (2D) and three-dimensional (3D) acoustic metasurfaces. At the same time, to further overcome the computational bottlenecks of traditional neural operators in large-scale grid divisions and high-dimensional data modeling, this paper focuses on predicting the scattering field only on the metasurface boundary. By combining this with the BEM to infer the full-field acoustic response, the degrees of freedom and grid size required for data generation are significantly reduced, thereby greatly improving the accuracy of both model training and inference. GRE-NO for introduces a GAN-based generative residual modeling to learn the full-field error distribution rather than deterministic point-wise residual regression, which outperforms conventional multilayer perception (MLP)-based residual correction networks in accuracy, convergence, and out-of-distribution generalization.

The main contributions of this work are summarized as follows:

1.    Generative residual enhancement of neural operators: We propose GRE-NO, a generative residual module built on top of the DeepONet + FNO framework, which explicitly learns the spatial distribution and statistical characteristics of the prediction error, improving prediction accuracy for boundary-based metasurface parameter analysis.

2.    Integration with BEM: By combining GRE-NO with BEM, our method predicts the acoustic scattering field only on the metasurface boundary, significantly reducing degrees of freedom and computational cost, while enabling efficient inference of the full-field response.

3.    High-accuracy cross-scale prediction under small data: GRE-NO effectively captures high-frequency and nonlinear error patterns, achieving accurate predictions even with limited training samples, which is critical for problems where generating large-scale datasets is computationally expensive.

4.    Comprehensive validation: We conduct extensive numerical experiments on both 2D and 3D metasurfaces, demonstrating that GRE-NO outperforms existing neural operator baselines, Including the standard FNO, DeepONet + FNO, and DeepONet + FNO + MLP Residual models.

The structure of this paper is as follows: Section 2 introduces the BEM formulation for acoustics and a detailed description of the GRE-NO method. Section 3 presents numerical examples on 2D and 3D metasurface parameter analysis, along with the corresponding results. Section 4 introduces and discusses the comparison between GRE-NO and other baseline neural network models. Section 6 concludes the paper and outlines directions for possible future work.

2  Method

To overcome the limitations in accuracy of some existing neural operators in cross-scale parameter analysis of metasurfaces, this paper proposes a new GRE-NO framework. Taking scattered acoustic field prediction as an example, GRE-NO introduces, for the first time, a residual generative modeling mechanism systematically into the neural operator architecture. Furthermore, by employing the BEM, the unbounded acoustic scattering problem is reduced to boundary-only prediction, thereby simplifying data complexity.

This section introduces the formulation of the acoustic BEM and the key features of the GRE-NO framework.

2.1 Boundary Integral Equation for the Acoustic Wave Problem

The scattering behavior of acoustic metasurfaces usually involves complex local resonance effects and far-field wave propagation, which is a typical unbounded domain acoustic scattering problem. Although traditional methods such as the Finite Element Method (FEM) and the Finite Difference Method (FDM) can provide accurate solutions for metasurface scattering problems, they require discretization of the entire computational domain and the use of absorbing boundary layers (e.g., perfectly matched layers, PML) to simulate unbounded conditions. This leads to significant increase of computational cost.

In contrast, BEM requires computation only on the structure boundary, naturally satisfying the unbounded radiation condition and significantly reducing the number of degrees of freedom. For example, in two-dimensional acoustic problems, the boundary integral method only requires discretization along the one-dimensional boundary, whereas in three-dimensional problems, it suffices to process the two-dimensional surface boundary. The fundamental idea of the BEM is to transform partial differential equations, such as the Helmholtz equation, into boundary integral equations.

Accordingly, our neural operator training is based exclusively on data defined on the metasurface boundary. To aid understanding, a brief summary of the BEM formulation is provided below.

In the frequency domain, linear acoustic waves are governed by the Helmholtz equation:

∇2ϕ+k2ϕ=0,∀x∈D(1)

where ϕ is the complex acoustic pressure, k is the wavenumber. On the boundary S, three types of boundary conditions are typically specified:

ϕ=ϕ¯; or q≡∂ϕ∂n=iωρvn; or ϕ=Zvn(2)

where ϕ¯ is the boundary sound pressure, q represents the normal derivative of the pressure, vn is the boundary sound velocity, Z is the acoustic impedance, i is the imaginary number, ω is angular frequency and ρ is density.

In addition, the radiation condition at infinity must be satisfied by the field:

limR→∞R(∂ϕ∂R−ikϕ)=0(3)

where R denotes the radius of an enclosing sphere around the scattering structure, and ϕ is the corresponding acoustic field.

By applying the Green’s second identity and the properties of the Green’s function, Eq. (1) can be reformulated into the Boundary Integral Equation (CBIE):

c(x)ϕ(x)=∫SG(x,y,ω)q(y)dSy−∫SF(x,y,ω)ϕ(y)dSy+ϕI(x)(4)

The Green’s function K(x,y,ω) depends on the spatial dimension of the problem:

G(x,y,ω)={i4H0(1)(kr),in2D14πreikr,in3D(5)

and

F(x,y,ω)=∂G∂n={−ik4H1(1)(kr)rjnjr,in2D(14πr2(ikr−1)rjnjr)eikr,in3D(6)

here, ϕI(x) is the incident field, Hn(1) denotes the Hankel function of the first kind, r=‖x−y‖, and nj is the unit normal at y.

The coefficient c(x) in Eq. (4) depends on the location of the point x, given as follows:

c(x)={1,x∈D12,x∈S(smooth boundary)(7)

In practical applications, the primary design objective of acoustic metasurfaces is to control the distribution of sound pressure ϕ(x) and the normal component of sound velocity q(x) on the scattering boundary. Once these quantities are specified on the boundary, the entire acoustic field, including both the near field and the far field, can be efficiently computed using the boundary integral equation. For a field point xf∉S, Eq. (4) yields:

ϕ(xf)=∫SG(xf,y,ω)q(y)dS(y)−∫SF(xf,yω)p(y)dS(y)+ϕI(xf)(8)

This equation shows that the complete acoustic field is uniquely determined once the boundary pressure and its normal derivative are known. However, in a conventional BEM framework, each change in geometric parameters or incident conditions requires reassembling and solving a new boundary integral equation. For acoustic metasurface design, especially in scenarios involving extensive geometric variations, broadband frequency analysis, or high-dimensional parameter spaces, this leads to thousands of repeated and computationally expensive BEM simulations. Such repetitive computations consume substantial time and computational resources, which significantly limits their applicability in inverse design, parametric studies, and large-scale optimization. Therefore, developing a model that can rapidly predict the boundary solution under varying parameters, and subsequently recover the full acoustic field through Eq. (8), is highly desirable for acoustic metasurface analysis.

2.2 Residual Compensation Mechanism of GANs

Against this background, neural operator methods such as the FNO and DeepONet have emerged as powerful tools for parameter-driven acoustic field modeling. These models offer end-to-end learning capabilities, strong generalization in high-dimensional functional spaces, and the ability to bypass repeated discretization and numerical solution procedures by learning the mapping from geometric parameters, material properties, and incident conditions to the corresponding acoustic fields. However, when applied to complex physical systems such as metasurfaces, neural operator models still face several challenges:

1.    Poor Out-of-Distribution (OOD) Generalization: The generalization performance of neural operators is highly dependent on the training data distribution. When new samples deviate from the training domain, the models may produce systematic errors even though they are capable of end-to-end predictions.

2.    Decreased Accuracy Near Boundaries and Singularities: Since neural operators often rely on global Fourier or kernel operations, significant prediction errors may occur near boundaries or in certain parameter regions (e.g., frequency transition points). Such deviations can lead to structural design failures or reduced accuracy in engineering applications.

3.    Difficulty in Quantifying and Correcting Errors: Neural operator predictions are effectively black-box outputs. Without a structured mechanism to characterize and compensate for residual errors on the boundary, even small deviations may accumulate when recovering the full acoustic field, limiting the robustness of small-sample or cross-scale predictions.

In view of the above limitations, we propose a residual modeling framework based on Generative Adversarial Networks to model and compensate for error distributions in complex physical systems or prediction models. The objective is to achieve accurate boundary predictions for acoustic metasurfaces using as few training samples as possible, thereby enabling rapid parameter analysis and full-field acoustic response computation via the BEM.

In this framework, two neural network models, the generator G and the discriminator D are trained simultaneously. The generator G aims to learn and simulate the true error distribution by sampling from a latent variable z∼pz and generating error samples G(z). The discriminator D, on the other hand, aims to distinguish whether the input error samples originate from the true error data distribution pdata or from the pseudo-error samples generated by G. The overall framework is illustrated in Fig. 1.

Figure 1: Schematic diagram of the generative residual learning framework for operator learning.

To avoid overfitting of the discriminator on a small-sample dataset, we alternately perform multiple discriminator updates and a single generator update in each training round. This strategy maintains the discriminator in a near-optimal state with respect to the current generated error distribution, while allowing the generator to gradually learn a more realistic error structure. This mechanism is analogous to maintaining an “error trajectory” during training, which guides the generator toward the target distribution even in the early stages, before full convergence. Ultimately, after training is complete, the generator can produce error samples that are difficult for the discriminator to distinguish, thus achieving realistic simulation and compensation of system residuals. Conventional residual correction for neural operators typically uses deterministic MLP networks to regress point-wise errors, which cannot capture spatial correlation or multi-scale error distribution. In contrast, our GAN-based generative residual module learns the statistical distribution of full-field residuals, enabling more robust compensation for systematic bias, especially in cross-scale and out-of-distribution cases.

The advantages of this method include:

1.    Compensation for accuracy blind spots: The generative residual mechanism can actively “sense” and focus on the error distribution of the neural operator, especially at boundaries or physical singularities. Through compensatory modeling, it significantly improves prediction accuracy and reliability, effectively enhancing global neural operator methods.

2.    Modeling flexibility: Unlike traditional discriminative networks, generative models are capable of capturing multimodal or non-Gaussian characteristics of residuals, making them more suitable for modeling complex error distributions and improving prediction fidelity.

3.    Strong architectural compatibility: The residual network can be seamlessly integrated into any primary neural operator architecture as an auxiliary module, without altering the original structure. Therefore, it can enhance the accuracy of any network, including physics-informed neural networks that typically struggle with precision.

4.    Improved robustness and generalization: By generatively reconstructing residuals and compensating the original outputs, the overall model’s adaptability to out-of-distribution data is improved. This enhances the stability of neural operators across various structures and frequency ranges.

2.3 GRE-NO Model Architecture

This section introduces the specific implementation method of GRE-NO. We retain the DeepONet architecture with FNO branches in the main model to preserve the structural advantages of neural operators for function-to-function mapping problems. Since our research focuses on predicting the scattered sound pressure in acoustic metasurfaces, which inherently involves wave propagation under complex boundary conditions and exhibits strong spatial frequency characteristics, it is both reasonable and necessary to employ FNO as the basic branch network. Therefore, FNO is used in both branches of the DeepONet architecture: one branch encodes the coordinate information (e.g., measurement point locations), and the other encodes functional information (e.g., the incident wave distribution and metasurface structural parameters). The final prediction of scattered sound pressure is obtained by fusing the outputs of both branches via an inner product.

In order to effectively capture the mapping relationships between complex functions, this study designs a dual-branch FNO-DeepONet structure based on the DeepONet framework. The basic idea of DeepONet is to learn in a general operator learning setting, where the input data consists of a set of function sensor values at discrete points, such as {y1,y2,…,ym}, and the goal is to learn the mapping from a function ΦI to the output p^s.

To handle structural and scale differences between the input function and spatial coordinates, DeepONet uses two network branches: one is the Branch Net, which encodes the input function u; the other is the Trunk Net, which encodes the spatial location y. The features extracted from the two branches are combined via an element-wise product to form the final output:

p^(y)=∑k=1nbk(ΦI)⋅tk(y)+b0(9)

where bk(ΦI) and tk(y) are the outputs of the branch and trunk networks, respectively, and b0 is a bias term. In this study, we further introduce a dual-branch network structure based on FNO within the DeepONet framework. Specifically, both the Branch and Trunk networks are composed of multi-layer FNO modules to fully extract multi-scale frequency domain features from the input functions and spatial coordinates. Each layer in the FNO is updated using the following formula:

vt+1(x):=σ(Wvt(x)+𝒦(a;φ)vt(x)),∀x∈D(10)

where 𝒦(a;φ) is a non-local convolution operator computed via the Fast Fourier Transform (FFT), defined as:

(𝒦(a;φ)vt)(x)=ℱ−1(ℱ(κφ)⋅ℱ(vt))(x)=ℱ−1(Rφ⋅ℱ(vt))(x)(11)

here, Rφ is a tensor containing the Fourier coefficients of the kernel function, which is learned adaptively through network training. This architecture effectively suppresses noise and extracts dominant physical features by preserving low-frequency components in the frequency domain, while also improving training and inference accuracy.

Therefore, the dual-branch FNO-DeepONet architecture combines the general operator approximation capabilities of DeepONet with the frequency domain modeling advantages of FNO. It is well-suited for modeling complex mappings in acoustic metasurface scattering problems and can capture the intricate relationships between high-dimensional structural parameters and scattered field responses. However, prediction accuracy may still be insufficient in some cases.

The loss function for the primary prediction task is defined as:

ℒmain=Eu,y∼𝒟‖p^y−py‖2(12)

Although the dual-branch FNO-based DeepONet framework can effectively learn complex multiscale operator mappings, it may still suffer from prediction bias under low-data training conditions due to limited data volume or insufficient coverage of high-dimensional features. To further enhance accuracy and robustness in such limited data scenarios, this work introduces a generative residual compensation mechanism to strengthen the modeling capability of neural operators.

Specifically, the output of the standard DeepONet is regarded as a first-order approximation of the true solution on the scattering boundary p(y), where the prediction error primarily arises from the insufficient capture of higher-order features. Therefore, we design an independent generative residual network that directly learns the residual between the primary prediction and the ground truth. During inference, this residual is used to compensate the original prediction, yielding a more accurate result.

Let p^(y) denote the prediction from the dual-branch FNO-DeepONet model and p(y) the ground truth. The residual is then defined as:

r(y)=p(y)−p^(y)(13)

The residual generator G takes as input the main network features, including the input y, and the main prediction p^(y), and then outputs a generated residual p^(y).

The adversarial training objective follows the original GAN formulation:

minGmaxDLadv=Er∼Pr[log⁡D(r)]+Er∼G[log⁡(1−D(r))](14)

where r^=G(y,p^(y)) is the generated residual. The final corrected prediction is given by:

pfinal(y)=p^(y)+r^(y)(15)

During training, the residual generator G is optimized by minimizing the mean squared error (MSE) between the true residual and the generated residual:

Lresidual=Eu,y||r(y)−G((y,p^(y)))||2(16)

To ensure good generalization of the residual generator without overfitting to noise, an appropriate L2 regularization term is added to the loss function to control the complexity of the generator’s parameters. Fig. 2 illustrates the framework of GRE-NO.

Figure 2: GRE-NO architecture diagram.

In addition, all of our data is sampled along the structure’s boundary. This not only conforms to the fundamental principle in physical field problems that boundary data determine the global solution, but also effectively reduces data dimensionality and storage costs. In the specific implementation, the GRE-NO model is used to predict the distribution of scattered sound pressure on the boundary of the metasurface structure, rather than regressing the full-field sound pressure throughout the computational domain.

Subsequently, the boundary scattered sound pressure predicted by GRE-NO is used as a boundary condition and input into the BEM. Taking advantage of BEM’s numerical strengths in solving wave problems, we can accurately compute the distribution of the entire sound field. Compared with traditional neural network approaches based on full-field supervision, our method significantly reduces the need for training samples and avoids the risk of overfitting that may arise when regressing full-field sound pressure in high-dimensional space.

Benefiting from this strategy of boundary prediction combined with reconstruction through mathematical and physical models, our aim is to achieve accurate prediction of broadband and complex acoustic fields using only moderate numbers of data samples. This is far fewer than the number of samples typically required by standard neural operators or deep networks.

3  Numerical Examples

To examine the effectiveness of the proposed GRE-NO method, a series of 2D and 3D scattering problems are used. In most problems, a plane wave with unit amplitude is incident on a set of rigid scatterers. Hard wall boundary conditions are imposed in all problems. In all the examples studied in this paper, the ground truth data is generated using the conventional boundary element method. In all numerical examples, train loss refers to the loss curve on the training data, test loss represents the loss curve on additional test data outside the training dataset.

3.1 Two-Dimensional Scattering Problem—Elliptical Scatterer

A hard-walled ellipse is subjected to a unit amplitude plane incident wave propagating in the positive x-direction. The geometry and related parameters are shown in Fig. 3, where “R1” denotes the semi-major axis of the ellipse and “R2” denotes the semi-minor axis. In this physical scenario, we explore two types of problems: frequency sweep problems with a fixed structure, and simulations of various structures at a fixed frequency. We first evaluate the prediction accuracy of GRE-NO for both problem types. Subsequently, we perform an ablation study to analyze the impact of different hyperparameters on the model performance.

Figure 3: Schematic diagram of an ellipse subjected to a plane wave of unit amplitude.

3.1.1 Prediction Accuracy Verification

The prediction accuracy is evaluated using two error metrics: Root Mean Square Error (RMSE) and relative Root Mean Square Error (rRMSE), which are defined as follows:

RMSE=1N∑i=1N(p^i−p¯i)2, rRMSE=1N∑i=1N(p^i−p¯ip¯i)2(17)

here, p¯i denotes the ground truth boundary pressure computed using the conventional BEM, p^i is the boundary pressure predicted by GRE-NO, and N is the total number of boundary points per sample. The rRMSE provides a more intuitive percentage-based interpretation of the error, while the RMSE is better suited for evaluating accuracy when the denominator is relatively small. All reported errors are presented in the rRMSE format. To ensure the statistical stability and repeatability of the model, all reported results in this study are obtained through five-fold cross-validation, and the error metrics are presented as “mean ± standard deviation”.

In the frequency sweep problem, R1 and R2 are kept constant. The frequency ranges considered are shown in Table 1. In the structural variation problem, the frequency is fixed, and the ranges of R1 and R2 are listed in Table 1 as well. All results in Table 1 are obtained using consistent GRE-NO parameters. Training samples are uniformly drawn from the full parameter space, while test samples are randomly selected from the remaining set. The Adam optimizer is employed with a learning rate of 0.0001, and the loss function is the L2 norm loss, defined as the mean squared difference between the predicted and true values. Based on the errors shown in Table 1, it is evident that GRE-NO performs well in both types of problems.

For the frequency-varying scenario, 1000 boundary solutions are generated using the conventional BEM. Of these, 700 are randomly selected for training and the remaining 300 for testing. Fig. 4 illustrates the acoustic pressure solution for a 2D elliptical scatterer predicted by GRE-NO. In this case, the analytical solution for the acoustic pressure is available and is used to compute the absolute prediction error:

ps=−p0∑n=0+∞εninJn′(ka)Hn(1)′(ka)Hn(1)(kr)cos⁡nθ−θd(18)

here, p0 is the pressure of the incident plane wave, r is the coordinate of the midpoint of the computational domain, and θ is the polar angle of the propagation direction of the incident wave, which is zero in this case, indicating the use of a first-order Hankel function. It should be noted that the training dataset does not include circular scatterers, but the alternative model can still provide good results, as shown in Fig. 5. Fig. 5a shows the loss curves of different datasets, where the training loss is the loss curve of the training dataset, the test loss is from a test dataset randomly selected from the training set, and test loss 0 is the test result of a dataset other than the training data. Fig. 5b shows the distribution of scattered sound pressure at the boundary points of a two-dimensional circular scatterer, where the horizontal axis represents the node index. The results show that the GRE-NO model can accurately predict the scattered sound pressure. Fig. 5c shows the prediction error of GRE-NO, and the comparison data is the analytical solution calculated by Formula (18).

Figure 4: Comparison of loss function curves for different hyperparameters of the GRE-NO model. (a) Loss function curves of different Fourier model components. (b) Loss function curves of different Fourier layers. (c) Loss curves of different generative residual network parameters.

Figure 5: Training results for a 2D circular scatterer. (a) Loss curve for the structural variation problem; (b) comparison of boundary pressures calculated by BEM and GRE-NO; (c) absolute error of the outer wavefield of the circular scatterer calculated by GRE-NO.

3.1.2 Ablation Studies

In this section, we examine the influence of various hyperparameters on the training performance of GRE-NO, including the Fourier mode number in the FNO branch, the number of Fourier layers, and the configuration of parameters within the GAN component. Two representative scenarios are considered: the frequency sweep problem and the structural variation problem. For the former, an elliptical geometry with a major semi-axis of 1.5, a minor semi-axis of 1, and a frequency range spanning from 1 to 4000 Hz is employed as the benchmark. For the latter, an elliptical geometry with a major semi-axis ranging from 0.5 to 3, a minor semi-axis from 0.2 to 2.7, and a fixed frequency of 500 Hz is utilized.

Fig. 4 illustrates the evolution of training loss under different hyperparameter configurations in the context of the structural variation problem. The dataset used in this case comprises 120 boundary grid points. Fig. 4a presents the results corresponding to three typical choices of Fourier mode numbers: 16, 32, and 64. It is observed that both the 16-mode (red) and 32-mode (blue) models exhibit rapid convergence and stable error reduction. Between these, the 16-mode configuration is preferable due to its lower computational cost. In contrast, when the number of modes increases to 64, the training loss initially decreases rapidly but plateaus after the second epoch, indicating potential overfitting. This behavior is likely to degrade generalization performance, as the model may learn non-generalizable or spurious patterns. Consequently, it is concluded that the number of Fourier modes should not exceed 32.

Fig. 4b examines the effect of varying the number of Fourier layers. It is evident that increasing the number of layers to 2 leads to a tendency toward overfitting. Therefore, a configuration consisting of a single Fourier layer and 16 modes is adopted in subsequent experiments. Fig. 4c investigates the sensitivity of the residual generator network to different parameter configurations (Parameter Sets 1, 2, and 3 vs. the Baseline). The baseline configuration, defined by a latent dimension of 64, a hidden layer dimension of 128, and three layers, is used as the reference. Three alternative configurations, as detailed in Table 2, are evaluated. The results indicate that all configurations eventually converge to comparably low training errors. This suggests that the residual generator is relatively robust to moderate changes in its architectural parameters, provided that they are set within a reasonable range.

While this subsection explores hyperparameter sensitivity, a structural ablation study comparing the base model and a deterministic residual correction (without adversarial training) is provided in Section 4. As detailed in later Fig. 6, the comparison between GRE-NO and the “DeepONet + FNO + MLP Residual” baseline shows explicitly the performance gain attributed to the generative adversarial mechanism.

Figure 6: Comparison of GRE-NO with other networks. (a) Relative error curve for the frequency scanning problem of phononic crystal scattering; (b) relative error curve for the structural change problem of phononic crystal scattering; (c) relative error curve for the frequency scanning problem of 3D spiral metasurface scattering; (d) relative error curve for the structural change problem of 3D spiral metasurface scattering.

3.2 Two-Dimensional Scattering Problem—Phononic Crystal Scatterer

Phononic crystals are a class of functional materials that employ artificially periodic structures to control the propagation of sound waves. They are widely used in applications such as bandgap engineering, acoustic waveguides, sound insulation, and sound energy manipulation. In this study, we focus on a two-dimensional phononic crystal unit cell with local resonance characteristics. The structure consists of an outer ring and a concentric inner ring, connected by a truncated narrow channel, forming a quasi-ring resonant configuration. This design incorporates not only the Bragg scattering effect arising from periodicity but also introduces a pronounced local resonance mechanism, which can generate a broadband acoustic bandgap within a specific frequency range.

When an incident acoustic wave excites the connecting channel of the unit cell, it induces local resonance within the inner ring, leading to energy localization and suppression of overall wave propagation. Near the resonant frequency, the structure exhibits strong wave-blocking capabilities. The resonance characteristics are primarily determined by the geometric parameters, with the size of the connecting channel having a decisive influence on the equivalent mass and stiffness.

3.2.1 Prediction Accuracy Verification

The geometry and associated parameters of the phononic crystal are illustrated in Fig. 7, where R1 and R2 denote the outer and inner radii, respectively, and h represents the width of the truncated channel. To assess the influence of structural variations on wave propagation characteristics, the values of R1, R2, and h are systematically varied to encompass a range of structural configurations.

Figure 7: Schematic diagram of a phononic crystal subjected to a unit amplitude plane wave.

Similar to the previous scenario, this study investigates both frequency-dependent and structure-dependent changes, with the objective of evaluating the modeling accuracy and generalization capability of the GRE-NO neural operator under varying input conditions. Within this framework, GRE-NO maps the spatial coordinates and incident frequencies of the boundary nodes to the corresponding acoustic pressure responses, thereby enabling prediction of the acoustic behavior of the entire Phononic Crystal (PnC) structure. As in the earlier example, the boundary solution is represented in the form of incident plane wave pressure, ensuring consistent input-output mapping across different cases.

For the frequency variation problem, a total of 1000 samples were generated with a sampling interval of 4 Hz, or the structural variation problem, 4150 samples were generated, each of which was analyzed using the traditional boundary element method. Among these, 70% of the boundary solutions were used for training, while the remaining 30% were randomly selected for testing. The training samples were uniformly distributed across the entire range of geometric parameters.

Due to the resonant characteristics of elastic waves, the boundary sound pressure near resonance exhibits rapid changes and differs significantly from the non-resonant pressure behavior. To accurately capture such complex relationships, it is generally necessary to sample more densely near the resonant frequencies. The resonant frequency of each elastic wave can be estimated by analyzing the pressure peaks across the frequency spectrum. Notably, GRE-NO is capable of maintaining high prediction accuracy without requiring specialized sampling around the resonant frequencies, demonstrating its superior generalization ability.

Table 3 presents the average values of RMSE and rRMSE for the test dataset, indicating that the training data provided strong overall accuracy. The results shown in Fig. 8a,b correspond to non-resonant and resonant PnC cases, respectively. A comparison of these results reveals that the prediction accuracy in the non-resonant case is comparable to that of the resonant case, indicating that no special sampling strategy is necessary.

Figure 8: Comparison of the prediction accuracy of GRE-NO at resonant and non-resonant frequencies in the acoustic scattering problem of phononic crystals.

Fig. 9 illustrates the prediction performance of the GRE-NO model for two types of problems: frequency sweep and structural variation. In the frequency sweep problem, the model demonstrates good convergence characteristics across different datasets. Fig. 9b further compares the relative errors of these datasets, all of which remain within 0.5%. Fig. 9c presents the prediction results at a frequency of 100 Hz, showing that the sound pressure values predicted by the model are highly consistent with those obtained using the boundary element method (BEM), thereby confirming the model’s reliability in handling frequency sweep problems.

Figure 9: Prediction results of GRE-NO in frequency sweeping problem and result change problem. (a) Comparison of training results for frequency sweeping problem; (b) comparison of relative errors of different datasets; (c) Frequency f = 100 Hz: Performance evaluation at a fixed frequency of 100 Hz; (d) comparison of training results for structural change problem. (e) Relative error analysis in structural change scenario; (f) changes of training and test loss of sample 3686 at 500 Hz with epochs, with the horizontal axis showing the loss distribution between node indices.

In the structural variation problem, the adaptive training results of GRE-NO across different topological configurations reveal that the predictions from various datasets converge effectively. The relative error analysis in Fig. 9e shows that the maximum relative error does not exceed 1%, reflecting the model’s strong generalization capability. Although the relative mean square error slightly increases for the structural variation case, this may be attributed to the increased complexity arising from cross-scale structural changes. Fig. 9f illustrates the prediction of sound pressure at the boundary nodes for sample 3686 at 500 Hz. The results indicate that the prediction error is slightly higher at the edge nodes compared to the central region, likely due to boundary effects; nonetheless, the error remains below 1%.

Fig. 10 compares the BEM solution, the GRE-NO prediction, and their error distribution for a phononic crystal scattering problem. After solving the boundary values, the sound pressure field is computed using Eq. (4). As shown in Fig. 10a, the BEM solution exhibits pronounced scattering effects near the slit channel of the phononic crystal, especially strong local amplification in the direction of the slit opening. Fig. 10b shows the sound pressure field predicted by the GRE-NO model, which is in high agreement with the BEM results, particularly in the primary scattering regions, where it accurately reproduces wave phenomena. Fig. 10c illustrates the error between the predicted and reference solutions, with the largest errors concentrated around the scatterer boundary, especially near the slit. This may due to the fact that these regions exhibit high pressure gradients and complex local resonance, which pose significant challenges for global operators. The generative residual module specifically targets these high-frequency discrepancies, and while it achieves high overall accuracy, the remaining small errors are inherently linked to these physically complex regions. However, the overall error magnitude remains below 0.05, indicating that GRE-NO achieves excellent near-field accuracy and captures high-frequency details effectively. These results demonstrate that the introduction of the residual generation module allows the model to learn localized structural discrepancies that are challenging for standard DeepONet + FNO frameworks, thereby improving predictive performance—especially for problems involving geometric singularities or strong local gradients.

Figure 10: Predicted and true sound pressure field with a phononic crystal scatterer under plane wave incidence. (a) True solution obtained by the boundary element method. (b) Sound pressure field predicted by GRE-NO. (c) Absolute error between the predicted and true sound pressure fields. The inset shows the geometry of the scatterer and the boundary value distribution.

3.2.2 Verification of the Generalization Ability of GRE-NO in Periodic Structures

Since phononic crystals are composed of periodically arranged units, even slight variations in their periodic parameters can significantly influence the resulting field distribution. To evaluate the generalization capability of the GRE-NO in periodic structures, this section investigates a representative periodic configuration of phononic crystals and assesses the model’s predictive accuracy under periodic conditions.

Similar to previous experiments, this study focuses on the scattered acoustic pressure response at the boundary of the periodic phononic crystal arrangement. The training dataset is generated using the conventional boundary element method and includes scattered field distributions corresponding to various periodic lengths. To assess the robustness of GRE-NO in extrapolation tasks involving periodic structures, the test data are deliberately selected from configurations that are not included in the training set.

The results are presented in Fig. 11. GRE-NO predictions within periodic structures exhibit convergence, and on the unseen test set, the relative error in scattered field prediction remains below 1.5%, indicating strong extrapolation ability and physical consistency.

Figure 11: The prediction results of acoustic pressure of phononic crystals with different period numbers. (a) The prediction loss curve of scattered acoustic pressure of phononic crystals with 1, 3 and 6 period numbers; (b) the prediction error curve of scattered acoustic pressure of phononic crystals with 1, 3 and 6 period numbers.

These findings demonstrate that GRE-NO maintains stable predictive performance under small-scale periodic perturbations. This generalization capability is primarily attributed to the residual enhancement mechanism, which enables adaptive modeling of localized field features. Moreover, the inherent global perception capability of the operator-based neural architecture helps prevent the overfitting commonly observed in purely data-driven models when applied to periodic structures.

3.3 Three-Dimensional Scattering Problem—Spiral Metasurface Scatterer

We observe that by longitudinally stretching the phononic crystal in the direction perpendicular to the plane and applying basic geometric transformations, it can be extended into a three-dimensional spiral metasurface. To further evaluate the accuracy and robustness of our model in more complex structural configurations, this section investigates three-dimensional chiral symmetric spiral absorbers. The proposed structure consists of a cylindrical body with a hollow core and a spiral path coiled around its surface, as illustrated in Fig. 12. This absorber is designed to exhibit anomalous sound absorption behavior and to suppress acoustic waves at specific target frequencies. Such capability is particularly important for applications involving frequency-selective noise control, such as active noise cancellation systems. The multi-parameter tunability of the spiral geometry enables a vast range of configuration combinations, each yielding distinct frequency responses. This versatility is essential for achieving advanced and precise control over broadband acoustic behavior, highlighting the potential of the proposed structure for next-generation metasurface-based acoustic manipulation.

Figure 12: Schematic diagram of the three-dimensional spiral metasurface and its related parameters.

3.3.1 Verification of Prediction Accuracy

The key design variables of a spiral metasurface include its pitch h1, width w1, and path thickness t1, as the acoustic absorption frequency is primarily influenced by these parameters. Therefore, in this study, we simulate wave scattering in spiral metasurfaces with varying lengths, widths, and thicknesses across a range of frequencies. The implementation of such complex geometries highlights the versatility and computational capability of the GRE-NO framework, demonstrating its superior performance in modeling intricate physical phenomena.

Fig. 12 illustrates the parameter definitions for a single spiral metasurface unit (SM). A total of 4000 samples with diverse structural parameters were generated, and the boundary acoustic pressure for each configuration was calculated using the traditional BEM. Of these, 3000 samples were allocated for training, while the remaining 1000 were randomly selected for testing. The training samples were uniformly distributed across the geometric parameter space to ensure comprehensive coverage.

This section does not include investigations involving assemblies of multiple three-dimensional metasurface units. Unlike phononic crystals, which rely on periodic arrangements to establish band structures, the functionality of metasurface designs is primarily governed by the resonance enhancement or attenuation behavior of individual unit cells at specific frequencies. The detailed parameter configurations and corresponding results are presented in Table 4.

Fig. 13 illustrates the prediction performance of the GRE-NO model in the acoustic scattering problem of spiral metasurface structures. As shown in Fig. 13a, for the 250th sample in the test set, the boundary sound pressure distribution predicted by the model closely matches the ground truth. The prediction error remains on the order of 10−5, indicating that the model effectively captures detailed variations in the sound field under complex geometries and achieves high prediction accuracy. Fig. 13b presents the training and testing loss curves, which exhibit a clear convergence trend. The training error decreases rapidly within the first 50 epochs and eventually stabilizes at an extremely low level. Simultaneously, the test losses stay close to the training loss, with no evident signs of overfitting, demonstrating the model’s strong generalization capability. Fig. 13c shows the evolution of the relative error during training. Both the training and test errors drop quickly and stabilize, with the final relative error maintained below 0.05. This further confirms the model’s robustness in high-accuracy scenarios.

Figure 13: Evaluation of the prediction performance of GRE-NO in the acoustic scattering problem of spiral metasurfaces. (a) Comparison of the true value of the boundary sound pressure of the 250th sample (left), the model prediction value (middle), and the prediction error (right); (b) the loss function curve of the model on the training set and the test set; (c) the change of relative error with the training epoch.

In summary, the GRE-NO model exhibits excellent modeling capability and numerical stability when handling three-dimensional spiral metasurface structures with high geometric complexity, offering a reliable approach for rapid prediction of complex sound fields. Figs. 14 and 15 display the results of the frequency sweep and structural variation problems, respectively. For each case, three representative samples are shown, including the BEM solutions, GRE-NO predictions, and the corresponding prediction errors. These visualizations intuitively demonstrate the high prediction accuracy of GRE-NO in complex three-dimensional structures.

Figure 14: Prediction results of spiral metasurface on a frequency scanning problem under complete data.

Figure 15: Prediction results of spiral metasurface on a structural parameter variation problem.

3.3.2 Missing Data Problem

This subsection investigates the effectiveness and robustness of the proposed GRE-NO model in handling incomplete training data, a common challenge in real-world acoustic measurement and simulation scenarios. The study focuses on predicting the boundary sound pressure in three-dimensional spiral metasurfaces, a representative multi-scale structure characterized by intricate geometry and strong spatial variability. Such complexity demands high accuracy and generalization capabilities from predictive models.

To rigorously evaluate model robustness, we constructed four distinct training datasets by randomly removing 0% (complete baseline), 5%, 10%, and 20% of the boundary sound pressure data points from the original high-resolution set. The model was trained on these incomplete datasets and evaluated against the complete ground truth.

Table 5 summarizes the prediction performance, measured by Relative Mean Squared Error, under different levels of data missingness. The results demonstrate remarkable robustness: with up to 10% of data randomly missing, the performance degradation is minimal, with rRMSE remaining on the order of 10−3. Even with a 20% missing rate, the error 10−2 remains within an acceptable range for many engineering applications. This defines a practical operational boundary for the model and provides a reliability benchmark for scenarios involving partial or corrupted experimental data.

Fig. 16 visually illustrates the data missing scenario and the model’s reconstruction capability. Subfigure (a) depicts the deliberately removed boundary data points, creating an incomplete input. Subfigure (b) shows the GRE-NO’s prediction based on this incomplete input, which is verified against the complete boundary data. The results confirm that GRE-NO can accurately reconstruct the full acoustic field even with significant portions of training data missing, highlighting its superior fault tolerance and ability to learn from sparse information.

Figure 16: GRE-NO prediction under incomplete training data. (a) Visualization of randomly removed boundary data points (missing proportion: 20%). (b) Predicted sound pressure field using the incomplete input from (a), compared against the complete ground truth.

3.3.3 Robustness Evaluation of Data Noise

To assess the generalization capability of the model under noisy conditions, we analyzed its performance by injecting Gaussian noise at varying intensities specifically, 0%, 1%, 5%, and 10% into the training data to simulate different levels of real-world data distortion. The corresponding prediction errors are summarized in Table 6.

The results indicate that as the noise level increases, the prediction error of the model exhibits a gradual, expected rise. With only 1% noise introduced, the model error 2.99e−4 remains within the same order of magnitude as the noise-free baseline 2.46e−4. This reflects the strong adaptability of GRE-NO to minor data perturbations. This characteristic underscore the model’s practical value in engineering applications.

When the noise level increases to 5%, the model error rises to 4.02e−3, indicating a noticeable decline in performance and suggesting that the noise intensity is approaching the effective robustness boundary of the model. Further increasing the noise to 10% leads to an error of 9.30e−3. Although the error progressively increases with higher noise levels, it remains below the order of e−2 even under 10% noise. This implies that the model does not fail due to data distortion, rather, it retains a reasonable level of physical predictive capability, highlighting its robustness in delivering meaningful predictions even when data quality is substantially compromised.

3.3.4 Three-Dimensional Periodically Arranged Spiral Metasurface

Building upon the investigation of individual spiral metasurfaces in the preceding section, we now extend our analysis to a more complex and technologically relevant configuration, i.e., the three-dimensional periodically arranged spiral metasurface. In practical acoustic engineering applications, metasurfaces are rarely employed as isolated units. Rather, they are typically arranged in periodic arrays to achieve enhanced wave manipulation capabilities over extended spatial domains. Such periodic arrangements introduce additional physical phenomena including Bragg scattering, bandgap formation, and collective resonance effects that cannot be observed in single-element configurations. The transition from a solitary metasurface to a periodic array substantially increases the computational complexity of traditional simulation methods. The unit cell interactions, long-range coupling, and multiple scattering events within the array create a highly coupled system where the acoustic response becomes a nonlinear function of both the individual element geometry and the array’s spatial arrangement. This complexity presents significant challenges for conventional BEM approaches, as the computational cost scales dramatically with the number of interacting boundary elements and the necessity to resolve interference patterns across multiple wavelengths. In this section, we demonstrate the capability of our GRE-NO framework to accurately model these intricate periodic systems while maintaining computational efficiency. We particularly focus on how the neural operator handles the additional degrees of freedom introduced by periodicity, captures inter-element coupling effects, and predicts array-scale wave phenomena that emerge from collective interactions.

In our investigation, each unit cell is characterized by three key geometric parameters (h1,t1,w1) that predominantly govern its acoustic performance. With each parameter discretized into merely two representative values, the complete combinatorial space for a six-unit periodic array encompasses (23)6=262,144 distinct configurations. Exhaustively evaluating this parameter space using conventional high-fidelity solvers, such as the Boundary Element Method, would demand prohibitive computational resources and time, rendering direct numerical exploration impractical for design optimization or sensitivity analysis.

In stark contrast, the GRE-NO framework demonstrates remarkable data efficiency in navigating this high-dimensional design space. By training on a dataset of only 16,000 samples, representing approximately 6% of the full combinatorial space. This efficient learning capability effectively decouples the cost of design evaluation from the exponential growth of the parameter space, transforming what would otherwise be an intractable computational burden into a manageable surrogate-based exploration. The specific results are shown in Table 7.

Fig. 17 illustrates the predictive performance of the proposed model for samples with structural variations. As shown in Fig. 17a, the predicted acoustic pressure field exhibits excellent agreement with the ground truth in both global distribution and local details. Under this magnified view, the model accurately reconstructs the acoustic field along the boundaries, while the remaining prediction errors are mainly localized in regions with high curvature or rapid spatial variation. The overall magnitude of the error remains small, demonstrating good generalization capability.

Figure 17: (a) Comparison of the true acoustic pressure field, the predicted field obtained by the neural solver, and the corresponding prediction error for sample 200. The scatter plots show the 3D spatial distribution of pressure values pt, where the discrepancy between prediction and ground truth is visualized through a color-coded error map; (b) training and testing loss curves for the structural change problem, showing the convergence behavior of the proposed model over 500 epochs on a logarithmic scale; (c) training and testing error curves, illustrating the reduction of global prediction error as training progresses. To facilitate visualization of the spatial variations within individual elements, the geometry is intentionally stretched along the y-axis; only the aspect ratio is modified, while all numerical values remain unchanged.

Fig. 17b presents the loss convergence for the structural variation task. The training loss decreases rapidly during the initial epochs and gradually stabilizes thereafter, while the testing loss remains of the same order as the training loss throughout the training process. This behavior indicates that the model exhibits stable optimization without noticeable overfitting. Fig. 17c further shows the evolution of the prediction error, where the testing error follows a similar decreasing trend as the training error, confirming that the model maintains strong predictive performance on unseen samples. Overall, the results demonstrate that the model retains high accuracy and stability even in scenarios involving complex structural variations.

4  Comparison of GRE-NO with Other Networks

To better evaluate the proposed GRE-NO framework, we clarify its relationship with several closely related state-of-the-art methods. The boundary-based FNO (Ref. [14]) achieves efficient parametric acoustic analysis using boundary-only input and Fourier neural operators, but does not employ residual enhancement or generative modeling. Residual-correction neural operators (Refs. [15,16]) improve accuracy via deterministic residual fitting, but they focus on point-wise error regression rather than full-field generative residual modeling. In contrast, GRE-NO innovatively integrates boundary-based operator learning with a GAN-based generative residual module to model the statistical distribution of prediction errors, leading to better accuracy, robustness, and out-of-distribution generalization.

In this section, we compare the performance of GRE-NO with several baseline network architectures. To ensure a fair comparison, all four methods are trained and tested on the same dataset, using identical batch sizes and optimization algorithms.

We first discuss the FNO, a neural network architecture based on the Fourier transform. The core idea of FNO is to learn mappings between functions in the frequency domain, enabling efficient prediction of partial differential equation solutions. Its architecture introduces Fourier layers to capture global spatial features and demonstrates strong generalization capabilities, especially suitable for physical modeling problems over continuous domains.

Next, we examine the DeepONet + FNO architecture, a unique network composed of two parallel FNO branches: the branch net, which processes functional inputs (such as physical parameters or initial conditions), and the trunk net, which handles spatial variables like coordinates. These two branches are fused via an inner product to predict the target function. DeepONet + FNO excels at learning functional mappings from input to output functions and is well-suited for modeling complex relationships in both low- and high-dimensional input spaces.

We also consider a standard residual compensation structure based on DeepONet + FNO, referred to as DeepONet + FNO + MLP Residual (DeepOnet + FNO + ResMLP). In this design, the residual compensation module is implemented using a multilayer perceptron (MLP), with skip connections that add the input directly to the output of the deep network. This allows the network to learn the residual mapping rather than directly fitting the target function. While the MLP-Residual shares conceptual similarity with GRE-NO, it implements a general residual correction mechanism.

Fig. 6 illustrates the training relative error curves for all networks on both the phononic crystal scattering problem and the 3D metasurface scattering problem. Results show that GRE-NO consistently outperforms the other three architectures under the same training data and number of iterations. GRE-NO exhibits the fastest convergence rate and the highest prediction accuracy. For the 2D frequency-sweep problem, all four neural networks achieve relative errors below 10%; however, GRE-NO reaches a relative error of 4% as early as the second iteration, highlighting the superior ability of generative models to accurately capture error distributions.

When addressing problems involving variations in two-dimensional structures, only GRE-NO achieves the desired prediction accuracy. The standard DeepONet performs the worst, with predictions nearly failing. The performance of both ResNet and FNO also degrades significantly. In contrast, GRE-NO maintains consistent performance on phononic crystal scattering problems, highlighting its ability to capture the underlying wave physics. This robustness may be attributed to the inclusion of the error probability distribution within the generative residual model.

GRE-NO enhances the primary operator model by introducing a generative residual network designed to learn the spatial distribution and statistical characteristics of prediction errors. Rather than directly approximating the target function, the residual network focuses on high-frequency perturbations and local discontinuities that the primary model cannot capture. This enables GRE-NO to establish a tighter link between physical consistency and numerical accuracy. Furthermore, by incorporating a generative adversarial mechanism, the residual distribution learning aligns more closely with the true error distribution, effectively mitigating overfitting or underfitting issues in highly nonlinear regions. As a result, GRE-NO improves prediction accuracy in cross-scale wave modeling tasks.

Notably, GRE-NO demonstrates exceptional predictive capabilities even with a very limited amount of training data. In this study, GRE-NO relies on only a few thousand training samples to accurately model scattering problems in phononic crystals with complex structural variations, achieving predictions comparable to high-fidelity simulation results. In contrast, under the same data scale, conventional network architectures such as FNO, ResNet, and DeepONet fail to produce reliable predictions. This finding clearly indicates GRE-NO’s significant advantage in small-sample learning, breaking through the reliance on large-scale data typically required by neural networks in complex physical systems. This leap in data efficiency not only reduces the computational cost of data generation but also provides a practical solution for physical modeling problems where data acquisition is limited, particularly in scenarios involving expensive simulations, large parameter spaces, and difficulties in discrete sampling.

In summary, each of the evaluated networks has its own strengths, but GRE-NO demonstrates superior fitting capacity and generalization performance by combining the benefits of generative residual learning and operator-based architectures. This makes it particularly suitable for acoustic parameter analysis tasks in high-dimensional parameter spaces.

5  Computational Efficiency

This section presents a quantitative comparison between the traditional BEM and the GRE-NO through a series of acoustic scattering cases, ranging from two-dimensional to three-dimensional geometries. The results reveal fundamental differences in computational paradigms and efficiency.

As shown in the accompanying Table 8, the computation time for BEM strictly adheres to the classical scaling laws inherent to its numerical formulation. The solution time increases from 1 s for a simple 2D case to 177 s for a complex 3D periodic case, representing an increase of nearly two orders of magnitude. This scaling is dictated by the O(N2) memory requirement and the O(N3) solution complexity associated with the dense system matrices generated by BEM. Consequently, for design optimization tasks requiring thousands of parameter evaluations, the total computational time can extend to several days (e.g., approximately 10 days for the 3D periodic case). This constitutes a fundamental bottleneck for the parametric exploration of metasurfaces. It is crucial to emphasize that this represents a limited-scale parameter study. As discussed in the previous section, for a periodic unit cell design with merely 3 design variables, each taking only 2 discrete values, the total number of combinations exceeds 260,000. In practical design workflows, relying on traditional numerical methods for exhaustive or even random search is computationally prohibitive. Designers are often forced to rely on intuition and trial-and-error, leading to inefficient processes with no guarantee of locating optimal solutions.

In contrast, the GRE-NO framework achieves a paradigm shift in efficiency by decoupling the computational load. Its “train once, infer many” strategy consolidates the primary cost into a one-time, offline training phase (e.g., approximately 108 min for the 3D periodic case). The online inference stage then consists of stable and highly efficient forward propagations. The key performance advantage lies in the fact that GRE-NO’s online inference time is completely decoupled from the discretization level (number of nodes) of the original physical problem, being determined solely by the fixed architecture of the neural network. For the most complex 3D periodic case, a single inference requires only 0.078 ms, representing an acceleration exceeding six orders of magnitude compared to a single BEM solve. When processing a batch of 4915 design samples, the total time for GRE-NO is merely 38.59 s, achieving an speedup of over four orders of magnitude compared to the approximately 869,955 s required by BEM. This capability makes real-time, interactive design exploration and large-scale parameter space scanning feasible.

An intriguing phenomenon observed in the analysis is that the per-sample prediction time for the 2D ellipse (0.0022 s) is slightly higher than that for the geometrically more complex 3D spiral (0.0014 s). We posit that this seemingly anomalous result stems primarily from the range of variation in the design parameters, rather than the geometric dimensionality of the problem. In the present study, the shape parameters for the 2D ellipse case vary over a meter scale with a relatively large range, leading to test samples that are more sparsely distributed and spaced farther apart in the normalized parameter space. Conversely, the dimension parameters for the 3D spiral case vary over a centimeter scale with a more concentrated range, causing test samples to more frequently reside within the high-probability density region covered by the training data. For a fixed, trained network, processing familiar inputs located in the core of the training distribution results in more direct and efficient activation patterns across the computational graph. In contrast, when input parameters lie at the edges of the training distribution or in wider-spaced regions, the network may need to engage deeper nonlinear expressive capacities for adaptation, microscopically increasing the complexity and execution time of the computational graph. Nevertheless, the efficiency fluctuation induced by parameter distribution is minimal and does not compromise the overall, absolute efficiency of GRE-NO’s millisecond-scale inference.

6  Conclusion

This paper introduces a new method for parametric acoustic wave analysis, GRE-NO, which has been validated through various numerical experiments by mapping problem boundaries and other physical parameters to the boundary solution space. These experiments include tests on acoustic wave scattering by rigid elliptical structures, phononic crystals, and spiral metasurfaces, covering both two-dimensional and three-dimensional cases. Notably, the required amount of training data remains on the order of thousands, even for complex three-dimensional metasurfaces. The results demonstrate that GRE-NO can be effectively applied to a wide range of structures and frequencies. Comparisons with other neural network-based models highlight its superior accuracy. Unlike previous studies, which primarily focused on modeling in finite domains, GRE-NO reduces the problem dimensionality by one and can model infinite-domain problems, further improving its data collection efficiency. Looking ahead, we plan to further integrate generative models with boundary element method, develop more physically informed neural architectures, and apply model reduction techniques to optimize the computational and memory efficiency for three-dimensional problems. While our current experiments focus on scattering problems, we are confident in the adaptability of GRE-NO and believe it can be extended to other complex problems, such as elastic waves and electromagnetic wave issues. In summary, GRE-NO is a promising framework for parametric wave analysis, particularly in infinite-domain settings involving complex geometries and frequency-sweeping tasks. Further research is required to fully explore its potential and address any remaining challenges. Despite the high accuracy achieved, this study primarily focuses on rigid boundary conditions. Future research could extend the GRE-NO framework to handle complex impedance boundaries and multi-physics coupling to further broaden its applicability in real-world acoustic environments.

Acknowledgement: This work was supported by the National Natural Science Foundation of China (Grant No. 12372198).

Funding Statement: This research was funded by National Natural Science Foundation of China (Grant number 12372198).

Author Contributions: Huilan Wu: Conceptualization; Methodology; Software; Data curation; Formal analysis; Visualization; Writing—original draft preparation. Yijun Liu: Conceptualization; BEM formulation, Funding acquisition; Supervision; Validation; Writing—review and editing. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author, upon reasonable request. The source code and the datasets generated for the GRE-NO framework during the current study are available in the https://github.com/hillllan.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Li J, Wang W, Xie Y, Popa BI, Cummer SA. A sound absorbing metasurface with coupled resonators. Appl Phys Lett. 2016;109(9):091908. doi:10.1063/1.4961671. [Google Scholar] [CrossRef]

2. Peng X, Ji J, Jing Y. Composite honeycomb metasurface panel for broadband sound absorption. J Acoust Soc Am. 2018;144(4):EL255–61. doi:10.1121/1.5055847. [Google Scholar] [PubMed] [CrossRef]

3. Wu H, Zhang H, Hao C. Reconfigurable spiral underwater sound-absorbing metasurfaces. Extreme Mech Lett. 2021;47:101361. doi:10.1016/j.eml.2021.101361. [Google Scholar] [CrossRef]

4. Cheng Y, Zhou C, Yuan BG, Wu DJ, Wei Q, Liu XJ. Ultra-sparse metasurface for high reflection of low-frequency sound based on artificial Mie resonances. Nat Mater. 2015;14(10):1013–9. doi:10.1038/nmat4393. [Google Scholar] [PubMed] [CrossRef]

5. Chen CT, Gu GX. Generative deep neural networks for inverse materials design using backpropagation and active learning. Adv Sci. 2020;7(5):1902607. doi:10.1002/advs.201902607. [Google Scholar] [PubMed] [CrossRef]

6. Rong J, Ye W. Multifunctional elastic metasurface design with topology optimization. Acta Mater. 2020;185(2):382–99. doi:10.1016/j.actamat.2019.12.017. [Google Scholar] [CrossRef]

7. Gu GX, Chen CT, Buehler MJ. De novo composite design based on machine learning algorithm. Extreme Mech Lett. 2018;18(36):19–28. doi:10.1016/j.eml.2017.10.001. [Google Scholar] [CrossRef]

8. Raccuglia P, Elbert KC, Adler PDF, Falk C, Wenny MB, Mollo A, et al. Machine-learning-assisted materials discovery using failed experiments. Nature. 2016;533(7601):73–6. doi:10.1038/nature17439. [Google Scholar] [PubMed] [CrossRef]

9. Bastawrous MV, Chen Z, Ogren AC, Daraio C, Rudin C, Brinson LC. A multiscale design method using interpretable machine learning for phononic materials with closely interacting scales. Comput Meth Appl Mech Eng. 2025;440(4):117833. doi:10.1016/j.cma.2025.117833. [Google Scholar] [CrossRef]

10. Li Y, Yin G, Yan G, Yao S. Forward-backstepping design of phononic crystals with anticipated band gap by data-driven method. Mech Syst Signal Process. 2025;224(4):111975. doi:10.1016/j.ymssp.2024.111975. [Google Scholar] [CrossRef]

11. Donda K, Zhu Y, Merkel A, Fan SW, Cao L, Wan S, et al. Ultrathin acoustic absorbing metasurface based on deep learning approach. Smart Mater Struct. 2021;30(8):085003. doi:10.1088/1361-665x/ac0675. [Google Scholar] [CrossRef]

12. Cheng H, Guo J, Zhang X, Ye W. Frequency-multiplexed transmitted-wave manipulation with multifunctional acoustic metasurfaces. Phys Rev Appl. 2023;20(3):034009. doi:10.1103/physrevapplied.20.034009. [Google Scholar] [CrossRef]

13. Kovachki N, Li Z, Liu B, Azizzadenesheli K, Bhattacharya K, Stuart A, et al. Neural operator: learning maps between function spaces with applications to PDEs. J Mach Learn Res. 2023;24(89):1–97. [Google Scholar]

14. Li R, Ye W, Liu Y. A boundary-based Fourier neural operator (B-FNO) method for efficient parametric acoustic wave analysis. Eng Comput. 2025;41(4):2393–410. doi:10.1007/s00366-024-02103-x. [Google Scholar] [CrossRef]

15. Cao L, O’Leary-Roseberry T, Jha PK, Oden JT, Ghattas O. Residual-based error correction for neural operator accelerated infinite-dimensional Bayesian inverse problems. J Comput Phys. 2023;486:112104. doi:10.1016/j.jcp.2023.112104. [Google Scholar] [CrossRef]

16. Jha PK. Residual-based error corrector operator to enhance accuracy and reliability of neural operator surrogates of nonlinear variational boundary-value problems. Comput Meth Appl Mech Eng. 2024;419:116595. doi:10.1016/j.cma.2023.116595. [Google Scholar] [CrossRef]

17. Lu L, Jin P, Pang G, Zhang Z, Karniadakis GE. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nat Mach Intell. 2021;3(3):218–29. doi:10.1038/s42256-021-00302-5. [Google Scholar] [CrossRef]

18. Haghighat E, bin Waheed U, Karniadakis G. En-DeepONet: an enrichment approach for enhancing the expressivity of neural operators with applications to seismology. Comput Meth Appl Mech Eng. 2024;420(1):116681. doi:10.1016/j.cma.2023.116681. [Google Scholar] [CrossRef]

19. Li Z, Kovachki N, Azizzadenesheli K, Liu B, Bhattacharya K, Stuart A, et al. Fourier neural operator for parametric partial differential equations. arXiv:2010.08895. 2020. [Google Scholar]

20. Lehmann F, Gatti F, Bertin M, Clouteau D. 3D elastic wave propagation with a factorized Fourier neural operator (F-FNO). Comput Meth Appl Mech Eng. 2024;420(1):116718. doi:10.1016/j.cma.2023.116718. [Google Scholar] [CrossRef]

21. Zhu M, Feng S, Lin Y, Lu L. Fourier-DeepONet: Fourier-enhanced deep operator networks for full waveform inversion with improved accuracy, generalizability, and robustness. Comput Meth Appl Mech Eng. 2023;416(6):116300. doi:10.1016/j.cma.2023.116300. [Google Scholar] [CrossRef]

22. Wagner JE, Burbulla S, de Benito Delgado M, Schmid JD. Neural operators as fast surrogate models for the transmission loss of parameterized sonic crystals. In: NeurIPS 2024 Workshop on Data-Driven and Differentiable Simulations, Surrogates, and Solvers; 2024 Dec 14–15; Vancouver, BC, Canada. [Google Scholar]

23. Kobayashi K, Daniell J, Alam SB. Improved generalization with deep neural operators for engineering systems: path towards digital twin. Eng Appl Artif Intell. 2024;131(1):107844. doi:10.1016/j.engappai.2024.107844. [Google Scholar] [CrossRef]

24. Baali H, Addouche M, Bouzerdoum A, Khelif A. Design of acoustic absorbing metasurfaces using a data-driven approach. Commun Mater. 2023;4(1):40. doi:10.1038/s43246-023-00369-0. [Google Scholar] [CrossRef]

25. Ryoo H, Jeon W. Broadband sound absorption using multiple hybrid resonances of acoustic metasurfaces. Int J Mech Sci. 2022;229:107508. doi:10.1016/j.ijmecsci.2022.107508. [Google Scholar] [CrossRef]

26. Zhao F, Zhang Z, Deng X, Feng J, Zhou H, Liu Z, et al. Atomic surface achieved through a novel cross-scale model from macroscale to nanoscale. Nanoscale. 2024;16(5):2318–36. doi:10.1039/d3nr05278h. [Google Scholar] [PubMed] [CrossRef]

27. Liu Y. Fast multipole boundary element method: theory and applications in engineering. Cambridge, UK: Cambridge University Press; 2009. [Google Scholar]

28. Liu Y, Rizzo FJ. A weakly singular form of the hypersingular boundary integral equation applied to 3-D acoustic wave problems. Comput Meth Appl Mech Eng. 1992;96(2):271–87. doi:10.1016/0045-7825(92)90136-8. [Google Scholar] [CrossRef]

29. Yang Y, Zhan Z, Liu Y. A novel damage identification algorithm by combining the boundary element method and a series connection neural network. Eng Appl Artif Intell. 2024;133(1851):108010. doi:10.1016/j.engappai.2024.108010. [Google Scholar] [CrossRef]

30. Honshuku Y, Isakari H. BEM-based fast frequency sweep for acoustic scattering by periodic slab. J Comput Phys. 2024;509(2):113046. doi:10.1016/j.jcp.2024.113046. [Google Scholar] [CrossRef]

31. Belibassakis K, Prospathopoulos J. A 3D-BEM for underwater propeller noise propagation in the ocean environment including hull scattering effects. Ocean Eng. 2023;286(1):115544. doi:10.1016/j.oceaneng.2023.115544. [Google Scholar] [CrossRef]

32. Kirkup S. The boundary element method in acoustics. Vol. 8. Preston, UK: University of Lancashire; 2007. [Google Scholar]

33. Gurbuz C, Kronowetter F, Dietz C, Eser M, Schmid J, Marburg S. Generative adversarial networks for the design of acoustic metamaterials. J Acoust Soc Am. 2021;149(2):1162–74. doi:10.1121/10.0003501. [Google Scholar] [PubMed] [CrossRef]

34. Brzin T, Jawed MK, Brojan M. Generative adversarial network-based inverse design of self-deploying soft kirigami composites for targeted shape transformation. Eng Appl Artif Intell. 2025;149(5):110417. doi:10.1016/j.engappai.2025.110417. [Google Scholar] [CrossRef]

35. Zhang H, Wang Y, Zhao H, Lu K, Yu D, Wen J. Accelerated topological design of metaporous materials of broadband sound absorption performance by generative adversarial networks. Mater Des. 2021;207(17):109855. doi:10.1016/j.matdes.2021.109855. [Google Scholar] [CrossRef]

36. Pan T, Pedrycz W, Yang J, Wang J. An improved generative adversarial network to oversample imbalanced datasets. Eng Appl Artif Intell. 2024;132(1):107934. doi:10.1016/j.engappai.2024.107934. [Google Scholar] [CrossRef]