BENO: Boundary-embedded Neural Operators for Elliptic PDEs

Q1: The paper only deal with Poisson equations with Dirichlet boundary conditions, which limits the effect on convincing readers that BENO is superior than baselines

A1: Thanks for your comment. We have generated new datasets on different boundary condition (i.e., Neumann BCs) under our existing setting with different corner cases. Typically boundary conditions can be categorized into Dirichlet and Neumann BCs, or the combination of the two. We believe that the addition of Neumann BCs to the manuscript makes the paper to be applicable to a much wider range of systems. We represent the supplementary experiment results as follows:

• Performances of our BENO and the compared baseline under Neumann boundary condition are shown below, which are trained on 900 4-corners samples and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix L.

• Performances of our BENO and the compared baseline under Neumann boundary condition are shown below, which are trained on 900 mixed samples (180 samples each from 5 datasets) and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix L.

Q2: How about more general elliptic equations?

A2: Thanks for the comment. We have generated new datasets on Darcy flow (mapping from coefficient to solution) under our existing settings. Performances of our BENO and the compared baseline on Darcy flow are shown below**,** which are trained on 900 4-corner samples (180 samples each from 5 datasets) and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix M.

Q1: It remains unclear to what extent the separation of the model into two terms is necessary without further empirical proof. Notably, both Branch 1 and Branch 2 share an identical architecture, differing only in their input, which comprises the forcing term and the boundary condition.

A1: Thanks for raising the question about different branches. Branch 1, with the boundary input set to zero, is posited to approximate the impact emanating from the interior, while Branch 2, nullifying the interior inputs, is conjectured to capture the boundary's influence on the interior. Specifically, we visualize the outputs of two branches which can fully validate our idea. The visualization is shown in the Appendix M in the updated manuscript, which indicates a discernible delineation of roles between the two branches.

Besides, you can also refer to BENO w/o. D in Table 6 Sec. 4.4 of the updated manuscript to see the effectiveness of our dual branch architecture

Q2: As illustrated in Figure 1, it seems that the encoding of boundary conditions is achieved using a shared transformer, while the forcing term employs two distinct MPNNs. Is this particular choice of architecture a necessary one?

A2: Thanks for highlighting the importance of examining the architecture of Transformers. In response, we have conducted an ablation study where we employed the MPNN architecture to encode boundary information, as an alternative to the Transformer model. The results of this study are presented below. They demonstrate that the MPNN architecture is less effective than the Transformer in encoding global boundary information, underscoring the superiority of the latter in this context.

Q3: Despite the acknowledged disadvantages of concatenating the boundary conditions alongside the other inputs, it is indeed surprising to see the catastrophic failure of recent advances, such as MP-PDE, as illustrated, especially in Figure 5. The authors have mentioned that the learning rate and the weight decay factor are fixed for all experiments through grid search, but it remains unclear with which model these hyperparameters are searched for. Furthermore, aligning optimization hyperparameters across all architectures may benefit certain models but harm the performance of others, and this aspect should be considered. Therefore, more justification should be issued to prove that other baselines do not work at all to support more concretely the claim of a performance gain.

A3: Thanks for the comment. In pursuit of a more equitable comparison, we have adopted additional strategies to enhance the baseline's performance. Firstly, we augmented the original MP-PDE with a virtual node linked to all boundaries (which has demonstrated to improve the performance, as is also supported by our experiments). Subsequently, we conducted an optimal hyper-parameter search to fine-tune the hyper-parameters. The ensuing results are presented below. Despite these efforts, a notable performance disparity remains when contrasted with our BENO approach. The baseline, even with these enhancements, does not successfully resolve the issue at hand.

Q4: From the visualization, it appears that there are visible artifacts aligned with boundaries. These artifacts may be introduced by including in the node embedding. This raises questions about the necessity of incorporating (dx, dy), as the Transform already appears to serve in learning a representation of the boundary shape

A4: Thanks for the comment. We have added ablation study that removes the dx, dy to validate its effectiveness, which is shown as follows. It is obvious that this feature is helpful for our PDE solving.

Q5: What motivates the use of different implementations for encoding the boundary condition (with Transformer) and the forcing term (with MPNNs)? Even if the goal is to learn a global representation for the boundary condition, couldn't MPNN achieve the same objective?

A5: Thanks for highlighting the importance of examining the architecture of Transformers. In response, we have conducted an ablation study where we employed the MPNN architecture to encode boundary information, as an alternative to the Transformer model. The results of this study are presented below. They demonstrate that the MPNN architecture is less effective than the Transformer in encoding global boundary information, underscoring the superiority of the latter in this context.

Q6: How are the 'distances to boundary ' defined? It appears to be more like an offset w.r.t the closest point on the boundary. If there are multiple closest points, how is the choice made?

A6: The calculation is based on the nearest point along the boundary. Given that the mesh points are arranged in a regular pattern, there is only one closest point for each calculation.

Dear Reviewer 113v,

We appreciate your time to review and the constructive comments to encourage our work. We want to leave a gentle reminder due to nearing the end of the discussion period.

We have tried to address all your concerns by providing more explanations and results. Please go over our response, and if you have additional questions, please let us know.

Thank you,

the Authors.

Q1: However, the paper's focus on elliptic PDEs seems somewhat restrictive, especially considering that it primarily deals with the simplest case of the Poisson equation with Dirichlet boundary conditions. This might give the impression that the novelty of the paper is somewhat limited.

A1: Despite the focus on elliptic PDEs, we believe this is an important and critical research that can be applicable to a wide range of natural phenomena and engineering applications.

• Importance of boundary conditions for any PDE systems. It is known that any PDEs can be characterized as parabolic, hyperbolic, and elliptic PDEs, which further can be categorized as initial-boundary value problem (IBVP), and initial value problem (IVP) and boundary value problem (BVP), respectively. For any of these PDEs, the treatment of boundary conditions plays a critical role in accurately modeling natural phenomena and engineering applications. While we focus on the application of BENO on elliptic PDEs, we believe that the boundary condition treatment and incorporation into the overall solver is applicable to other PDE systems.
• Applicability of elliptic PDEs to the real-world phenomena. The utility of elliptic PDEs is not limited to steady problems, e.g., we can obtain unsteady solution of nonlinear PDEs by repeated solution of elliptic PDEs. For instance, the pressure-Poisson equation (which is an elliptic PDE) for incompressible flows can be applied to unsteady phenomena. Various plasma simulations uses Poisson solver that is time dependent. Helmholtz equation is used in a variety of phenomena, including electromagnetic radiation, seismology, and acoustics. These systems can have complex geometries, e.g., airfoil, propulsion devices, plasma sources for semiconductor manufacturing, with different types of materials (e.g., dielectric, conducting, catalytic, chemically reactive) and inhomogeneous boundary values (e.g., dependent on space and time). Thus, we believe that demonstration of BENO on elliptic PDEs has immense impact to a large number of applications. This introduction is updated to emphasize this point. Besides, we have revised the manuscript to delete the statements that elliptic PDEs are only applicable for steady-state phenomena.

Furthermore, to validate the generalization ability of our BENO, we also supplement more experiments on other scenarios in the following part.

Besides, to validate the generalization ability of our BENO, we also supplement more experiments on other scenarios in the following part.

Q7: In the paper, it is mentioned that other methods overfit in challenging boundary conditions, but a more explicit explanation of what this means would be helpful. Are there graphs or figures that can demonstrate this concept?

A7: Thank you for bringing attention to the issue regarding the failure of other methods. We have updated our manuscript to reflect a more precise terminology, changing 'overfitting' to 'extremely hard to learn'. This revision enhances the accuracy of our statements. The evidence for this adjustment is evident from Figure 4 in Section 4.2 and Figure 5 in Appendix I.3. These figures clearly illustrate that the baseline methods struggle to converge the loss during both training and testing phases. Additionally, there is a markedly high prediction error when compared to the ground truth, further underscoring the challenges these methods face.

Q8: It would be beneficial to provide a clearer explanation of what the red and orange parts in Figure 1 represent. What does the red part signify, and what about the orange part? How do they differ?

A8: The red and orange hues in Figure 1 are chosen from the color-bar for our graphical representation, which effectively illustrates the boundary values. The redder the area, the higher the boundary value it represents, whereas the more orange the area, the lower the boundary value. We have added this information in the caption of Figure 1 in the updated manuscript.

Q9: Regarding Table 1, as mentioned in the Strengths section, is it impossible for FNO to handle 5.? Is it also impossible for GKN to handle 5.? Is there no prior research in the field of operator learning models that can handle 5.?

A9: To the best of our knowledge, no work has stepped in solving elliptic PDEs with complex boundary shape and inhomogeneous boundary value together. As for all of the baselines, we utilize the same information as input for fair comparison and it is obvious that existing vanilla baselines cannot deal with the complex setting, which just demonstrates that a new architecture design is needed.

Q10: In the Problem Setup (2. Problem Setup), it is not clear whether BENO can be applied to 3D or higher-dimensional domains. While the use of graph neural networks suggests that it might be possible, further explanation would be helpful.

A10: We have refined our setup in Section 2 to enable its generalization to more diverse scenarios.

Q11: In the context of 3.1 Motivation, where branch1 is interpreted as the interior and branch2 as the exterior from a Green function perspective, it's important to know whether the output of Branch1 precisely reaches 0 at the boundary and whether the output of Branch2 reaches exactly "g" at the boundary. Is this a process of approximating "g," or is the model capable of obtaining the exact boundary values of "g"?

A11: Thanks for raising the question about different branches. Branch 1, with the boundary input set to zero, is posited to approximate the impact emanating from the interior, while Branch 2, nullifying the interior inputs, is conjectured to capture the boundary's influence on the interior. Specifically, we visualize the outputs of two branches which can fully validate our idea. The visualization is shown in the Appendix N, which indicates a discernible delineation of roles between the two branches.

Q12: In Section 3.2, a more detailed explanation of the definition of would be appreciated. Furthermore, the part around this section involving attributes is somewhat challenging to comprehend in written form. It might be a good idea to exclude this part if it ultimately relates to the explanation in Section 3.3.2.

A12: We have refined the manuscript in alignment with your feedback.

Q13: Are dx, dy, and dc in Section 3.3.1 uniquely determined for each node? Also, why is dc necessary?

A13: Thanks for the comment. dx, and dy are uniquely determined for each node as the node feature. We have added ablation study that removes the dx, dy to validate its effectiveness, which is shown as follows. It is obvious that this feature is helpful for our PDE solving.

Q4: Furthermore, the experiments in complex domains appear to be limited to similar, complex domains, such as the "5 different corner elliptic dataset." It would be interesting to explore whether the BENO method can be extended and applied to different scenarios, such as domains with circular holes in the center of a square or corners at non-right angles (or more simple, circle domain as in (Lotzsch et al., 2022)). By conducting experiments in more diverse domains and presenting results in a variety of scenarios, the paper could better showcase the novelty and versatility of the BENO method.

A4: Thanks for the positive feedback. Previous work i.e., GNN-PDE, utilizes methodology on circular holes but may have neglected the treatment of complex physical boundary conditions (including combination of Dirichlet and Neumann boundaries, which can be inhomogeneous), which we believe is the key innovativeness of our work. Generalization of the method (e.g., circular, spherical, unstructured meshes) is reserved for future work.

Q5: Therefore, it would be beneficial to generalize the limitations of BENO, as described in this paper. For instance, in terms of boundary conditions, there are relevant discussions in [1], and graph-based methods have addressed complex domains and a variety of equations and boundary conditions, as seen in [2] and [3], which leveraged Finite Element Method (FEM). Referencing these works might provide valuable insights.

A5: Thanks for pointing out these references. We have added these references in the revised version of 5. Related Work section for a boarder audience.

Q6: Furthermore, there is some curiosity regarding whether the results of the baseline methods were adequately compared. It appears from Appendix F that the internal and boundary grids for the baseline methods are distinguished using one-hot encoding. It's worth considering if this method provides the fairest basis for comparison. As mentioned in the paper, FNO suggests the use of Fourier continuation for different boundary conditions and domains ('Non-uniform and Non-periodic Geometry' part in [4]). Therefore, it raises the question of whether the paper's approach is indeed the best for addressing complex scenarios. Additionally, there is a model for operator learning called DeepONet [5], which can be applied directly to challenging boundary conditions and complex domains. It would be interesting to know whether the paper conducted a comparison with this model.

A6: Thank you for your comment. In pursuit of a more equitable comparison, we have adopted additional strategies to enhance the baseline's performance. Firstly, we augmented the original MP-PDE with a virtual node linked to all boundaries. Subsequently, we conducted an optimal hyper-parameter search to fine-tune the hyper-parameters. The ensuing results are presented below. Despite these efforts, a notable performance disparity remains when contrasted with our BENO approach. The baseline, even with these enhancements, does not successfully resolve the issue at hand.

And it should be noted that the boundary conditions utilized in the related work are all homogeneous. We have also added the DeepONet as the baseline, and implement experiments under our main setting, which is trained on 900 4-corners samples and tested on 5 datasets. The results are shown as below and it is obvious that our BENO surpass the baseline by a large margin.

Q2: For example, can the proposed BENO method be applied to Poisson equations with Neumann boundary conditions or mixed Robin boundary conditions? The current model design appears to be tailored to the Dirichlet boundary conditions, which might be considered quite restrictive.

A2: Thanks for your comment. We have generated new datasets on different boundary condition (i.e., Neumann BCs) under our existing setting with different corner cases. Typically boundary conditions can be categorized into Dirichlet and Neumann BCs, or the combination of the two. We believe that the addition of Neumann BCs to the manuscript makes the paper to be applicable to a much wider range of systems. We represent the supplementary experiment results as follows:

• Performances of our BENO and the compared baseline under Neumann boundary condition are shown below, which are trained on 900 4-corners samples and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix L.

• Performances of our BENO and the compared baseline under Neumann boundary condition are shown below, which are trained on 900 mixed samples (180 samples each from 5 datasets) and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix L.

Q3: In this context, the paper experimentally investigates a single type of operator, the one between (f, g) and the solution u. However, it's intriguing to consider whether other types of solution operators could also be learned (for instance, as seen in the FNO paper with Darcy flow, where the coefficient could be an input to the operator, or the initial condition of PDE can be an input).

A3: Thanks for the comment. We have generated new datasets on Darcy flow (mapping from coefficient to solution) under our existing settings with different corner casess. Performances of our BENO and the compared baseline on Darcy flow are shown below**,** which are trained on 900 4-corner samples (180 samples each from 5 datasets) and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix M.

Q14: On page 5, is "B" ultimately equal to "B^1"? It would be beneficial to provide a more detailed explanation of where "t" changes and what its range is.

A14: t is the step idx ranging from 1 to 5 in our experiments. And B^t represents the boundary embedding in the t-th message passing step.

Q15: The definition of MAE for Table 2 is present in the appendix but not in the main text.

A15: Thanks for the catch. In the updated manuscript, we have added the definition of MAE for Table 2 to the main text.

Q16: In Section 4.3.2, the explanation of training at 32x32 and testing at 64x64, is this capability attributed to BENO's use of graph neural networks, or are there other features enabling this super-resolution? A more detailed explanation would be appreciated.

A16: This capability is primarily attributed to the GNN architecture used in BENO. The GNN's capability to handle super-resolution stems from its approximation of the Green's function. The Green's function describes how the value at one location influences the solution at the current position. This characteristic makes GNN particularly suitable for super-resolution tasks, as it adeptly models the interplay and impact of various spatial data points on each other.

Q17: In Appendix K, when looking at Figure 6, it appears that the results for "Ours" are different from other solution profiles, forming square blocks with discontinuities. What is the reason or cause behind this phenomenon?

A17: The key factor is the narrow range of numerical values in shown data. In such a scenario, even minor fluctuations in the predicted values near the corners can result in the values falling outside of a given scale interval on the color-bar. This creates the visual effect of discontinuous blocks, as opposed to a smoother gradient typically seen with a wider range of values.

Dear Reviewer MqG8,

We appreciate your time to review and the constructive comments to encourage our work. We want to leave a gentle reminder due to nearing the end of the discussion period.

We have tried to address all your concerns by providing more explanations and results. Please go over our response, and if you have additional questions, please let us know.

Thank you,

the Authors.

Q1: The grid sizes for the domain are quite small (32x32, 64x64) when compared with FNO. I would like to see if the approach generalizes to larger and more practical grid sizes as well. This is especially important given the sequence length constraints that come with using transformers.

A1: Thanks for the comment. We have conducted additional super-resolution experiment, where the models are trained on 64x64 resolution, and are tested on irregular domains at a high 128x128 resolution, pushing the boundaries of our BENO's applicability. The results are shown as below. Specifically, BENO consistently shows superior performance by a large margin.

We see that our model still outperforms the MP-PDE baseline by a large margin, at this high resolution.

Q2: It would also be useful to get some additional information about training and inference times for BENO and baselines.

A2: Thank you for your comment. The average time taken for BENO and the graph-based baseline methods during a single epoch, computed over five iterations, is listed below. We can see that the modest increase in computational time of BENO is justified by the substantial performance improvements offered by the integration of the Transformer and dual-branch GNN structures. The added complexity of these components contributes to the enhanced capability of the BENO model, making it a feasible and efficient choice in scenarios where advanced performance is crucial.

Q6: Can the BENO model be used to generate physically meaningful insights and interpretability?

A6: The innovativeness of BENO is in accounting for the boundary condition treatments in the PDE solutions obtained from boundary-embedded GNN. Handling of the boundary conditions is a requirement for any PDE equations. Elliptic PDEs are good test problems to test such boundary condition effects (cf. boundary value problem). In addition, without accurately capturing the boundary conditions, the solutions for the PDEs obtained from a solver are not physically meaningful and cannot be applied to real systems. Hence, we believe that BENO’s strengths are in producing an acceptable tolerance with much faster computational speed, which generate physically meaningful results and enhance the interpretability of the solutions.

Q7: There are some latest related operator works that should be referred to for a boarder audience:

A7: We have added [Xiao, Xiongye, et al. "Coupled Multiwavelet Neural Operator Learning for Coupled Partial Differential Equations."], [Gupta, Gaurav, et al. "Multiwavelet-based operator learning for differential equations."], and [Gupta, Gaurav, et al. "Non-linear operator approximations for initial value problems."] in the revised version of Related Work for a boarder audience.

Q3: How can the dataset used to evaluate the BENO model be expanded to include more diverse scenarios and improve the generalization performance of the model?

A3: Thanks for your comment. We have generated new datasets on different boundary condition (i.e., Neumann BCs) type and different types of solution operators (i.e., Darcy flow) under our existing setting with different corner cases. Typically boundary conditions can be categorized into Dirichlet and Neumann BCs, or the combination of the two. We believe that the addition of Neumann BCs to the manuscript makes the paper to be applicable to a much wider range of systems. We represent the supplementary experiment results as follows:

• Performances of our BENO and the compared baseline under Neumann boundary condition are shown below, which are trained on 900 4-corners samples and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix L. We see that our method significantly outperforms the MP-PDE baseline in this setting.

• Performances of our BENO and the compared baseline under Neumann boundary condition are shown below, which are trained on 900 mixed samples (180 samples each from 5 datasets) and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix L. Again, we see that our method significantly outperforms the MP-PDE baseline in this setting

• Performances of our BENO and the compared baseline on Darcy flow are shown below**,** which are trained on 900 4-corner samples and tested on 5 datasets under relative L2 norm and MAE separately. The whole comparison table with all the baselines is shown in Appendix M.

Q4: Would the authors be able to provide the number of parameters to quantify each model's complexity?

A4: We would like to present the parameter count for our BENO model. Specifically, the Transformer component, dedicated to encoding the boundary, comprises 0.8M parameters, while each branch of the GNN-based network encompasses 1.14M parameters. This indicates that the complexity of our approach is entirely manageable and within an acceptable range.

Q5: How can the BENO model be adapted to handle time-dependent PDEs or other types of dynamic systems?

A5: As we responded in A2, even some elliptic PDEs are designed to handle time-dependent solutions and the importance of boundary condition treatment remains same for other dynamic systems. We have modified our introduction to highlight this point.

Q1: Lack of comparison with more SOTA works like multiwavelet neural operator.

A1: Thanks for your insightful feedback. We have added the Multiwavelet-based Operator as the baseline, and implement experiments under our main setting, which is trained on 900 4-corners samples and tested on 5 datasets. The results are shown as below and it is obvious that our BENO surpass the baseline by a large margin.

Q2: Focusing on solving elliptic PDEs with complex geometries and inhomogeneous boundary values may limit the applications.

A2: Despite the focus on elliptic PDEs, we believe this is an important and critical research that can be applicable to a wide range of natural phenomena and engineering applications.

• Importance of boundary conditions for any PDE systems. It is known that any PDEs can be characterized as parabolic, hyperbolic, and elliptic PDEs, which further can be categorized as initial-boundary value problem (IBVP), and initial value problem (IVP) and boundary value problem (BVP), respectively. For any of these PDEs, the treatment of boundary conditions plays a critical role in accurately modeling natural phenomena and engineering applications. While we focus on the application of BENO on elliptic PDEs, we believe that the boundary condition treatment and incorporation into the overall solver is applicable to other PDE systems.
• Applicability of elliptic PDEs to the real-world phenomena. The utility of elliptic PDEs is not limited to steady problems, e.g., we can obtain unsteady solution of nonlinear PDEs by repeated solution of elliptic PDEs. For instance, the pressure-Poisson equation (which is an elliptic PDE) for incompressible flows can be applied to unsteady phenomena. Various plasma simulations uses Poisson solver that is time dependent. Helmholtz equation is used in a variety of phenomena, including electromagnetic radiation, seismology, and acoustics. These systems can have complex geometries, e.g., airfoil, propulsion devices, plasma sources for semiconductor manufacturing, with different types of materials (e.g., dielectric, conducting, catalytic, chemically reactive) and inhomogeneous boundary values (e.g., dependent on space and time). Thus, we believe that demonstration of BENO on elliptic PDEs has immense impact to a large number of applications. This introduction is updated to emphasize this point. Besides, we have revised the manuscript to delete the statements that elliptic PDEs are only applicable for steady-state phenomena.

Furthermore, to validate the generalization ability of our BENO, we also supplement more experiments on other scenarios in the following part.