HyPINO: Multi-Physics Neural Operators via HyperPINNs and the Method of Manufactured Solutions

We thank the reviewer for taking the time to carefully evaluate our work and for the feedback. We appreciate your insights and address your concerns below.

[Q1:] While our approach builds on prior work, the contributions of HyPINO go significantly beyond a simple combination thereof. Neural operators such as FNOs and CNOs typically require large datasets of high-fidelity simulation data for training. The introduction of the physics-informed loss in PINO reduced this dependency significantly. However, some supervised data is still necessary, as training solely with physics-informed objectives is often unstable due to spectral bias and optimization challenges [1]. This remaining dependency on simulation data is problematic when designing multi-physics solvers, as generating and combining data for diverse physical problems in a single dataset is challenging. This is why we propose to leverage the method of manufactured solutions (MMS) to obtain an (potentially) infinite number of diverse training samples. To the best of our knowledge, HyPINO is the first neural operator trained on a dataset spanning such a diverse set of linear, 2D PDEs with mixed Dirichlet and Neumann boundary conditions on complex domains.

This infinite dataset also allows us to revisit the previously underexplored HyperPINN design, which is a more data-intensive neural operator than CNOs or FNOs because it generates a continuous target function parametrized by a PINN as opposed to the grid-based outputs of CNOs and FNOs. The benefits of HyperPINNs are that (i) the output PINN is continuous and can be evaluated at arbitrary points on the domain, (ii) exact analytical derivatives can be obtained directly via automatic differentiation, and (iii) the computation of these derivatives is computationally inexpensive because the gradients for computing the residual only need to be backpropagated through the lightweight PINN rather than the large Swin Transformer. While previous work has recognized these benefits and used HyperPINNs as backbones of their neural operators, their scope was typically restricted to minor variations of a single equation or physical system (e.g. different source terms for the same PDE, or the same equation on multiple geometries). Through the combination with the diverse, infinite MMS-dataset we are able to unlock the potential and scale HyperPINNs to true multi-physics neural operators.

Furthermore, while HyPINO shows strong zero-shot performance for a diverse set of benchmark PDEs, consistently outperforming other baseline operators, it also naturally supports task-specific fine-tuning. This is done by taking the weights generated by HyPINO as initialization and training the resulting PINN with standard physics-informed objectives using analytical derivatives to compute and minimize the residual. This task-specific fine-tuning is impossible with most other neural operator designs because they do not produce a parameterized solution that is separate from the model.

Finally, we introduce the iterative refinement process, a training-free method that leverages superposition to improve the predictions of physics-informed neural operators for linear PDEs. Given an initial PINN solution from HyPINO, we can compute its PDE residual (i.e. how much it violates the PDE). This residual function can be fed back into the hypernetwork to generate a “delta-PINN” representing a correction to the original solution. Adding the delta-PINN’s predictions to the original PINN’s predictions yields an updated solution that cancels out much of the previous residual. This procedure can be iterated multiple times, thus creating an ensemble of self-correcting networks. This refinement does not require any retraining of the hypernetwork. It is a post-processing technique that uses the trained neural operator model repeatedly in inference mode. To the best of our knowledge, this iterative residual refinement idea with delta-PINNs is also entirely novel and possibly even applicable to other physics-informed neural operator designs, such as PINO, to improve their predictions on linear PDEs.

[1] Goswami, Somdatta, et al. "Physics-informed deep neural operator networks." Machine learning in modeling and simulation: methods and applications. Cham: Springer International Publishing, 2023. 219-254.

[Q2:] We agree that second-order optimizers such as L-BFGS or Gauss-Newton methods were shown to improve convergence in physics-informed training. In fact, many state-of-the-art PINN setups adopt a two-stage scheme with Adam followed by a second-order optimizer. However, in our experiments, the goal was not to reach state-of-the-art performance on each benchmark PDE, but to isolate and assess the impact of HyPINO’s initialization quality. To do so, we deliberately chose a minimalistic and widely used training configuration so that improvements could be directly attributed to the initialization. As shown in Figure 4, HyPINO-initialized PINNs (a) start with lower error, (b) have well-behaved convergence curves, and (c) reach lower errors faster than random and Reptile initializations. While we recognize that techniques such as Sobolev losses, dynamic loss weighting, adaptive sampling, and second-order optimizers can improve PINN performance, our aim was to evaluate HyPINO initializations in a clean and controlled setting. Including such enhancements would have shifted the focus from the quality of the initialization itself and introduced confounding factors that would have made a direct comparison difficult. That said, our method is fully compatible with more advanced PINN training strategies. If desired, these can easily be added on top during downstream adaptation. While we did not run experiments with second-order optimizers or other advanced PINN techniques, the benefits of HyPINO initializations are likely to carry over, as these strategies are complementary and largely independent of the initialization scheme.

[Q3:] Our training objective indeed combines both physics-informed and supervised terms when analytical solutions are available via MMS. Specifically, the loss used for training is $J=\lambda_R J_R+\lambda_D J_D+\lambda_N J_N + J_S$, where $J_R$ is the residual loss over interior collocation points, $J_D$ and $J_N$ are Dirichlet and Neumann boundary losses, and $J_S$ is a Sobolev supervised objective computing a Huber loss on function values, first, and second derivatives between the generated PINN and the analytical solution at sampled collocation points. This is detailed in Section 3.4 (Eq. 8). The Sobolev loss [2] was shown to significantly improve the accuracy and generalization of neural networks, including PINNs [3]. However, computing higher-order derivatives required for the Sobolev loss is typically computationally expensive and thus impractical in many applications. In our setting, the necessary derivatives are already obtained as part of the residual evaluation for the physics-informed loss, which allows us to incorporate the Sobolev loss without additional computational overhead.

[2] Czarnecki, Wojciech M., et al. "Sobolev training for neural networks." Advances in neural information processing systems 30 (2017).

[3] Son, Hwijae, et al. "Sobolev training for physics informed neural networks." arXiv preprint arXiv:2101.08932 (2021).

We thank the reviewer again for their comments and hope to have addressed their concerns satisfactorily. If so, we kindly request the reviewer to upgrade their assessment accordingly, while being at their disposal for any further questions about our paper.

[Q1, P1] We thank the reviewer for raising this point and would like to clarify that all baseline models were trained or fine-tuned on our MMS-generated dataset. As detailed in Section 4.2:

• U-Net, which uses the same encoder as HyPINO but replaces the hypernetwork with a convolutional decoder, was trained from scratch using MMS-labeled data and supervised loss.
• Poseidon, a pretrained neural operator, was fine-tuned on MMS-labeled data, also with a purely supervised objective.
• PINO was trained from scratch using the same hybrid objective (supervised + physics-informed) and curriculum as HyPINO.

Model performances are reported in Table 1 (Section 4.3). We hope this clarifies that the evaluation of established neural operators on our MMS-generated data is not only included, but a core part of our experimental design.

We would also like to comment on the "(classic) method of manufactured solutions" remark, which implicitly conveys that the MMS-based part of our work is not novel or not a significant contribution. While it is true that MMS is a well-established technique for verifying the correct implementation of numerical solvers and their convergence behavior, our contribution is to show that MMS can be used not just to validate, but for training neural operators capable of generalizing across a broad range of 2D PDEs, including linear elliptic, parabolic, and hyperbolic equations with mixed boundary conditions and complex geometries while also being able to quantify the solution error. This required a substantial effort in building a scalable and diverse data-generation pipeline that produces realistic samples and ensures that the resulting neural operator is not only effective on synthetic data, but also capable of transferring to real-world PDE benchmarks, as demonstrated in our evaluations. The data-generation pipeline will be made public so that other researchers can generate datasets for training and evaluating their physics-informed neural operators.

[Q1, P2] We agree that applying the residual-based refinement strategy to other physics-informed neural network designs is an interesting exercise. Therefore, we ran an additional experiment where we applied the same refinement strategy to PINO. The table below shows the performance across benchmarks, where superscripts indicate the number of refinement steps:

We observe consistent improvements with refinement, especially on HZ, PS-C, and PS-G. The drop in performance on PS-L is similar to what we saw with HyPINO and is likely related to the small target values. We plan to explore these findings further in future work.

Overall, these results support the idea that residual-based refinement is applicable to physics-informed neural operators in general and not limited to our specific architecture. If accepted, we would be happy to include this table and a short discussion in the appendix.

We also thank the reviewer for pointing out the connection to multi-stage networks. While, as mentioned by the reviewer, those methods train each stage independently, both approaches use iterative, residual-based improvements. We will briefly discuss this line of work and add a citation in Section 3.5.

We thank the reviewer for taking the time to carefully evaluate our work and for the constructive feedback. We appreciate your insights and address your points below.

[Q1:] We agree that resolution-equivariance is a key property for neural operators. While the output of HyPINO is a continuous PINN that can be evaluated at arbitrary spatial coordinates, the input PDE parametrization (source function and boundary masks/values) is discretized on a fixed-size grid (224×224) to match the Swin Transformer’s input resolution. Following prior work [1], this limitation can be mitigated by demonstrating test-time resolution invariance when varying the input grid resolution and resizing it to (224x224). We performed this ablation on the Helmholtz benchmark (HZ) by changing the source function resolution between 28 and 448:

Between resolutions of 56 and 448, SMAPE varied by less than 0.3, which shows approximate invariance. Only at very coarse resolutions (28x28) does the performance start deteriorating. This experiment can be added to the appendix.

[1] Herde, Maximilian, et al. "Poseidon: Efficient foundation models for pdes." Advances in Neural Information Processing Systems 37 (2024): 72525-72624.

[Q2:] In HyPINO, Swin blocks are used to encode the discretized inputs (i.e., the source fields, domain geometry, and boundary conditions) into a latent representation. To condition this representation on the PDE operator, we apply FiLM layers after each Swin block. These layers modulate the intermediate latent features conditioned on the PDE coefficients $c \in \mathbb{R}^5$. This is a lightweight and common mechanism for global conditioning, consistent with prior work in neural operators [1], and better suited than cross-attention in our setting, because the conditioning signal is a 1D vector rather than a sequence of tokens.

In the decoder, we adopt a hypernetwork-style design: pooled latent features are passed through MLPs to generate the weights and biases of the target PINN.

The Swin encoder, the FiLM-conditioning, and the hypernetwork decoder serve complementary roles, which makes isolated ablation of individual components less straightforward. However, while designing the architecture, we experimented with different Transformer backbones, MLP widths, and activation functions. We found these architectural choices to have relatively minor influence on overall performance compared to the quality of the dataset and the choice of loss functions, such as the addition of the physics-informed loss, whose effectiveness is evidenced by the superior performance of HyPINO over the UNet baseline that shares the same encoder. We therefore mostly adopted standard architectural settings from prior work.

[Q3:] We would like to emphasize that HyPINO does not backpropagate through any external numerical PDE solver. Our framework is entirely neural: HyPINO is a hypernetwork that maps a PDE specification $(L, f, g, h)$ to the weights $\theta^\star$ of a small target PINN $u_{\theta^\star}$. This PINN approximates the PDE solution as a continuous function over the domain. Residuals are computed by applying the differential operator $L$ to $u_{\theta^\star}$ via automatic differentiation, and the training objective is to minimize the residuals.

At inference time, our residual-based refinement strategy also avoids solver-level differentiation. We compute the residual of the current PINN and pass it through the hypernetwork in forward mode only to obtain a corrective delta-PINN. The updated prediction is then formed by summing the PINN outputs. Only the small PINNs are differentiated to compute residuals whereas the hypernetwork remains frozen.

Regarding computational costs, to obtain the loss of a predicted PINN $u_{\theta^\star}$, we evaluate it on a set of collocation points and compute the necessary first- and second-order derivatives to obtain the residuals. These higher-order derivatives are computed efficiently via automatic differentiation of the small PINN only (3 layers, width 32). The resulting residual loss is then backpropagated to update the hypernetwork. Importantly, this gradient only flows through the Swin-based encoder at the final step after residuals have been computed. Thus, while residuals require evaluating derivatives up to second order, the computationally expensive hypernetwork is only differentiated once per batch.

[Q4:] Thank you for pointing out these relevant works. We will include and discuss them in the related works section of the final version.

We thank the reviewer for taking the time to carefully evaluate our work and for the constructive feedback. We appreciate your insights and address your points below.

[W2:] All baselines were fine-tuned on our synthetic dataset using a similar training budget (same number of steps, batch size, and optimizer settings). While the UNet and Poseidon are trained without the physics-informed objective, PINO was trained with the same settings as HyPINO. We agree that the FNO architecture in PINO was never intended for the type of PDE parametrization used in this work (i.e., binary indicator grids for the boundaries), which might make it less effective than a Swin Transformer. Conversely, HyPINO must predict a full set of PINN weights, which is more challenging than predicting a solution grid. We acknowledge the difficulty of direct comparison in the evaluation Section 4.3 on line 250.

[W3:] It is correct that the iterative refinement process relies on linearity. That said, some non-linear PDEs can be linearized, making the method also applicable to those.

[Q1:] Our framework can be extended to 3D, non-linear, and higher-order PDEs with only minor modifications, such as adding non-linearities like $u u_x$ to the set of selectable terms when sampling the differential operator, or adding a third variable and changing the dimensionality of the inputs to the Swin Transformer and target PINNs. However, more computational resources would be necessary to train such a model. We want to highlight that, to the best of our knowledge, HyPINO is the first neural operator trained on such a diverse set of linear, 2D PDEs with mixed Dirichlet and Neumann boundary conditions on complex domains.

While our current framework focuses on linear PDEs, this still captures a broad and impactful class of PDE problems in engineering and the natural sciences, often sufficient for many modelling applications or as first steps into non-linear modelling approaches. Moreover, we have experimented with linearization techniques (e.g. the Cole–Hopf transformation for the 1D Burgers equation) as a practical extension strategy, which enables HyPINO to handle certain nonlinear PDEs. However, we agree this work should trigger novel research directions, such as including non-linear examples directly in the training pipeline.

[Q2:] Indeed, the part about constructing the binary mask for the boundary conditions could use some additional information. We sample continuous boundary points directly on the CSG-defined shapes and need to convert them into a grid-based discretization for the Swin Transformer. For each point on the boundary, we find the enclosing grid cell it lies in and set the four corner points of that cell to 1 in the binary mask (indicating that the boundary passes through this region). For example, a boundary point at $(0.5, 0.5)$ activates the grid points $(0,0)$, $(1,0)$, $(0,1)$, and $(1,1)$. This gives the binary masks $M_g$ and $M_h$. The corresponding boundary values are stored in $V_g$ and $V_h$ at the same locations. This can be clarified in the final version.

[Q3:] We agree that discretization-invariance is a key property for neural operators. While the output of HyPINO is a continuous PINN that can be evaluated at arbitrary spatial coordinates, the input PDE parametrization (source function and boundary masks/values) is discretized on a fixed-size grid (224×224) to match the Swin Transformer’s input resolution. Following prior work [1], this limitation can be mitigated by demonstrating test-time resolution invariance when varying the input grid resolution and resizing it to (224x224). We performed this ablation on the Helmholtz benchmark (HZ) by changing the source function resolution between 28 and 448:

Between resolutions of 56 and 448, SMAPE varied by less than 0.3, which shows approximate invariance. Only at very coarse resolutions (28x28) does the performance start deteriorating. This experiment can be added to the appendix.

[1] Herde, Maximilian, et al. "Poseidon: Efficient foundation models for pdes." Advances in Neural Information Processing Systems 37 (2024): 72525-72624.

I am overall satisfied with the rebuttal by the authors, but I am disappointed that the reference below, which came to my attention after I had written my review, was not cited or mentioned in the related works section:

Jae Yong Lee, SungWoong Cho, Hyung Ju Hwang, "HyperDeepONet: learning operator with complex target function space using the limited resources via hypernetwork", ICLR 2023.

I think the authors should include this reference and position their contribution in reference to that. I would like to keep the score as is.

We thank the reviewer for taking the time to carefully evaluate our work and for the constructive feedback. We appreciate your insights and address your points below.

[W1:] While it is true that most of the experiments focus on steady-state PDEs, our evaluation also includes the 1D wave equation (WV), a time-dependent hyperbolic PDE. HyPINO demonstrates competitive zero-shot performance on this benchmark (MSE: 2.9e-1), outperforming UNet and Poseidon, and matching PINO.

[W2:] The reviewer is correct in stating that the current work is limited to PDEs in (at-most) two space dimensions. However, our framework is very general and can be readily extended to 3-D PDEs, which we aim to consider in future work.

[W3:] Thank you for highlighting these relevant works. We will include and discuss them in the related works section of the final version.

[Q1:] Indeed, the choice and sampling procedure of collocation points used to evaluate the generated PINNs could be explained further. The target PINNs are evaluated on 4,096 randomly sampled points, where 2,048 lie within the domain, and 2,048 are sampled on the boundaries. The reasons for this are twofold: first, sampling on the entire grid (over domain and boundary) would yield 224*224 collocation points, which would not fit into memory. Second, by sampling explicitly on the boundary, we ensure better coverage and a stronger training signal for Dirichlet and Neumann boundary conditions. If we would only rely on points from the sampled grid, very few would lie exactly on the boundary, especially for complex shapes, which could lead to weak or noisy training signals for the boundary conditions. This sampling process will be clarified and more detailed in a camera-ready version of the paper, if accepted.