Calibrated Physics-Informed Uncertainty Quantification
Authors: Vignesh Gopakumar
Ander Gray
Lorenzo Zanisi
Timothy Nunn
Daniel Giles
Matt Kusner
Stanislas Pamela
Marc Peter Deisenroth
Venue: ICML 2025 poster (2025)
Keywords: Surrogate Models
Uncertainty Quantification
Neural-PDE
Physics-Informed
Conformal Prediction
PDE Residuals
Nuclear Fusion
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: Calibrated uncertainty quantification of neural-PDE solvers using physics residual errors as a non-conformity score for conformal prediction.

Abstract: Simulating complex physical systems is crucial for understanding and predicting phenomena across diverse fields, such as fluid dynamics and heat transfer, as well as plasma physics and structural mechanics. Traditional approaches rely on solving partial differential equations (PDEs) using numerical methods, which are computationally expensive and often prohibitively slow for real-time applications or large-scale simulations. Neural PDEs have emerged as efficient alternatives to these costly numerical solvers, offering significant computational speed-ups. However, their lack of robust uncertainty quantification (UQ) limits deployment in critical applications. We introduce a model-agnostic, physics-informed conformal prediction (CP) framework that provides guaranteed uncertainty estimates without requiring labelled data. By utilising a physics-based approach, we can quantify and calibrate the model's inconsistencies with the physics rather than the uncertainty arising from the data. Our approach utilises convolutional layers as finite-difference stencils and leverages physics residual errors as nonconformity scores, enabling data-free UQ with marginal and joint coverage guarantees across prediction domains for a range of complex PDEs. We further validate the efficacy of our method on neural PDE models for plasma modelling and shot design in fusion reactors.

Calibrated Physics-Informed Uncertainty Quantification
Vignesh Gopakumar' 12
Ander Gray? Lorenzo Zanisi² Timothy Nunn Daniel Giles'
Matt J. Kusner
Stanislas Pamela" Marc Peter Deisenroth
Abstract
Numerical Solver
FALSE
Output
Simulating complex physical systems is cru- cial for understanding and predicting phenomena across diverse fields, such as fluid dynamics and heat transfer, as well as plasma physics and struc- tural mechanics. Traditional approaches rely on solving partial differential equations (PDEs) using numerical methods, which are computationally ex- pensive and often prohibitively slow for real-time applications large-scale simulations. Neural PDEs have emerged efficient alternatives to these costly numerical solvers, offering significant computational speed-ups. However, their lack of robust uncertainty quantification (UQ) limits de- ployment in critical applications. We introduce a -ggnosic,,, -nysscs-informee conformal pre- diction (CP) framework that provides guaranteed uncertainty estimates without requiring labelled data. By utilising physics- based approach, we can quantify and calibrate the model's inconsisten- cies with the physies rather than the uncertainty arising from the data. Our approach utilises con- volutional layers finite-difference stencils and leverages physics residual errors as nonconfor- mity scores, enabling data-free UQ with marginal and joint coverage guarantees across prediction domains for range of complex PDEs. We further validate the efficacy of our method on neural PDE models for plasma modelling and shot design in fusion reactors.
IC, BC, Geometry, Physics
Training Data
Surrogate Model
UQ
Figure I. Neural PDE framework: Neural PDE solvers use data from traditional numerical solvers to quickly approximate PDEs across various conditions (shown by black arrows). To ensure reli- ability, these models incorporate uncertainty quantification (UQ) methods. If the predicted етгоr exceeds a coverage threshold c, numerical solver is utilised, further adding the training data; otherwise, predictions
used output (shown by red arrows), This paper focuses developing new UQ method assessing confidence in neural PDE models, emphasising the grey-shaded region of the framework.
celli et al., 2024), but their computational demands and environmental impact pose significant challenges (Horwitz, 2024). Machine learning approaches have emerged as effi- cient alternatives for modelling/emulating PDEs (Bertone et al., 2019; Karniadakis al., 2021), successfully deployed across weather forecasting (Kochkov al., 2024; Meyer et al., 2022: Giles al., 2024), fluid dynamics (Jiang et al., 2020; Pfaff et al., 2021), and nuclear fusion applications (Poels et al., 2023; Carey et al., 2024; Gopakumar & Samad- dar, 2020). Neural PDE solvers provide rapid approxima- tions, but present a critical cost-accuracy trade-off. While generating outputs consistently, their solutions may violate physical constraints or produce misleading results with high confidence (Gopakumar et al., 2023a). A typical neural- PDE framework (see Figure 1) trains surrogate models on numerical simulation data to predict the evolution of PDE under various conditions, with uncertainty quantification (UQ) methods reverting to numerical solvers when predic- tions fail coverage thresholds. However, current UQ meth- ods lack statistical guarantees (Zou et al., 2024), require extensive simulation data (Gopakumar et al., 2024a), require architectural modifications (Abdar et al., 2021).
1. Introduction
Numerical PDE solvers are essential tools in scientific and engineering simulations (Danabasoglu et al., 2020; Giudi-
'Centre for Artificial Intelligence, University College Lon- don *Computing Division, UK Atomic Energy Authority "Heudiasyc Laboratory 'Polytechnique Montréal Milaa Que- bee Al Institute. Correspondence to: Vignesh Gopakumar <v.gopakumar @ @uck.ac.uk>.
To address these limitations, we propose a framework com- bining PDE residuals over neural PDEs with Conformal Pre- diction (CP) to provide uncertainty estimates that guarantee coverage. Our approach evaluates Physics Residual Errors
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Calibrated Physics-Informee UQ
(PRE) from neural PDE solver predictions and performs calibration using marginal and joint CP formulations. Our method provides statistically valid coverage within the resid- ual space (Vovk et al., 2005), offering enror bounds against the violation of physical conservation laws. The framework is model-agnostic (applicable to all types of surrogate mod- els) and does not require additional data generated from expensive numerical simulations. It yields interpretable uncertainty bounds indicating the model's physics-guider inconsistencies, addressing the (over)confidence issue of neural PDE solvers (Zou et al., 2024). Our contributions:
2019) and has found application in optimising neural op- erators (Li et al., 2024) and soft/hard enforcement of the physical constraints to deep learning models (Du et al., 2024; Chalapathi et al., 2024). However, they have rarely been used as a tool for providing UQ to surrogate models, and where they have found application, UQ remained uncali- brated (Zhu et al., 2019). The majority of the literature in UQ for neural PDE solvers has been looking at Bayesian methods, such as dropout, Bayesian neural networks, and Monte Carlo methods (Geneva & Zabaras, 2020; Zou et al., 2024; Psaros et al., 2023), which lack guarantees or are computationally expensive.
Calibrated physics-informed UQ: A novel physics- informed nonconformity metric using PDE residuals. This quantifies uncertainty through physical inconsistency rather than training data variation, providing input-independent prediction sets while relaxing exchangeability restrictions.
3. Background
3.1. Neural PDE Solvers
Consider the generic formulation of a PDE modelling the spatio-temporal evolution of n field variables u E R across range of initial conditions:
Marginal and Joint CP: Our approach guarantees cover- age bounds both marginally (univariate per dimension) and jointly (multivariate across the entire prediction domain), en- abling the identification of volatile predictions and creating a peteecioo-acceptanee criteria.
== Di(u)= + Dx(u) = 0, XeΩ, tE [0,
(1)
w(X, = 9 XE
(2)
u (X,0) - X).
(3)
2. Related Work
Here, X defines the spatial domain bounded by Q, [0, temporal domain, and Dx and D; the composite operators of the associated spatial and temporal derivatives. The PDE is further defined by the boundary condition y and initial condition a, which can be parameterised by X. The set of solutions of field variables is expressed as u E U.
T] the
Recently, CP, as method of performing UQ, has been gain- ing popularity for usage with spatio-temporal data (Sun, 2022). Several works have explored the inductive CP frame- work for spatial and sequential data (Stankeviciute et al., 2021; Xu & Xie, 2021; Xu et al., 2023), including in the operator space (Ma et al., 2024). Gopakumar al. (2024a) extends the marginal-CP framework to pre-trained as well as to fine-tuned surrogate models for physical system mod- elling across an chinite-dimensiong setting. Alternatively, error bounds for PDE surrogates have been devised by Gray et al. (2025) using set propagation to project the singular value decomposition of the prediction error onto the predic- tion space.
Neural PDE solvers as surrogate models learn the behaviour governed by Equation (1) using a parameterised neural net- work NW Starting from the initial conditions, the network is trained to solve the spatio-temporal evolution of the fields given by Ω U [o, T). Neural operators NOp are a special class of neural networks that learn the operator mapping from the function space of the PDE initial conditions a E A to the function space of solutions u E U. A neural operator for solving an initial-value problem can be expressed as
Marginal
Prediction
Inaak Cenatioo
Plytis Reskalal
u = =Oo((AA), - -oo (M)
Neural POE
Physica
Residual
Prediction
Conformal
X LE [0,
(4)
Figure 2. Schematic of physics-informed uncertainty quantifica- tion workflow, Initial conditions generate neural PDE predictions autoregressively, over which physics residual errors are estimated. Calibration via marginal and joint conformal prediction yields er- ror bars pointwise for marginal-CP CP and domain-wiee for joint-CP.
A Fourier Neural Operator (FNO) is an autoregressive neural operator that learns the spatio-temporal evolution of PDE solutions by leveraging the Fourier transform the ker- nel integrator (Li et al., 2021). The field evolution is learned using tunable weight matrices of the network, parameterised directly in the Fourier space of the PDE solutions.
The usage of PDE residuals under the guise of Physics- Informed Machine Learning (PIML) (Karniadakis et al., 2021) was made popular as an optimisation strategy for Physics-Informn Informed Neural Networks (PINNs) (Raissi al.,
Since CP and our extension of it provide a post-hoc measure of quantifying the uncertainty of a neural PDE, remains
2

---

Calibrated Physics-Informee UQ
agnostic to model choice and training conditions. Consider- ing the model independence of our approach, we restrict our experiments to modelling PDEs with FNO. The FNO is chosen due to its cost-accuracy trade-off and efficiency as demonstrated by de Hoop et al. (2022) and Gopakumar et al. (2023b). CP over a range of neural-PDE solvers has been applied by Gopakumar et al. (2024a), who also demonstrate that the coverage guarantees are upheld irrespective of the model choice, further necessitating us to not experiment with various model architectures.
Table 1. Overview of nonconformity metries AER, STD, and PRE and their corresponding score functions and prediction sets.
Score Function (8) Prediction Sets (C*)
AER
(1F(X))
xi))
f(Xn+1) 2
(1ioor " (o(x)) 1-1
STD
fu(Xnti) # (X
(ID((XNN))
PRE
3.2. Conformal Prediction
Conformal prediction (CP) (Vovk et al., 2005; Shafer & Vovk, 2008) is a statistical framework that addresses the accuracy of a predictive model. Consider a machine learning model f: X' + y trained on a dataset (X,, Y)) that can be used to predict the next true label Ya+1 at query point X CP extends the point prediction P Y to a
restricts their application to domains where sufficient data exists a priori or can be easily obtained
4. Physics Residual Error (PRE)
prediction set Co, ensuring that
We introduce a novel data-free nonconformity score based directly on the PDE for surrogate models. The Physics Residual Error (PRE) defined as the PDE residual (Saad Schultz, 1986) estimated over the discretised PDE solution obtained from the surrogate model. For an abstract PDE as in Equation (1), the PDE residual is the evaluation of the composite differential operator D over a field(s) u under the influence of an external force b
P(Y.+1 E CC)21-00
(5)
This coverage guarantee, function of the user-defined con- fidence level o, holds irrespective of the chosen model and training dataset. The only condition is that the calibration samples and the prediction samples are exchangeable. Tradi- tional inductive CP partitions the labelled data into training and calibration sets (Papadopoulos, 2008). The performance of the model on the latter, measured using a Ronconformity score, is used to calibrate the model and obtain prediction sets.
D(u) b= 0,
The PDE residual is treated as a score function by taking its LI norm as indicated in Table 1. While well-defined PDEs have solutions obeying Equations (1) to (3), numerical so- lutions often fail to converge the true solution (Pinder, 2018). Neural PDEs, trained on approximate numerical data, are further prone to non-convergent predictions. In numerical analysis, the norm of the PDE residual is often used as a criterion for stability, convergence, and accuracy (Iserles, 2009). The PRE typically represents the violation of conservation laws associated with the physical system. Using the residual error as nonconformity score quantifies the neural PDE solver's non-convergence to the physical ground truth of the PDE. By further using conformal pre- diction, we obtain coverage bounds over this conservative, residual space with guaranteed coverage.
Conventionally, nonconformity scores act the model pre- dictions and a labelled dataset (Kato et al., 2023). For deterministic models, they are often formulated as the Ab- solute Error Residual (AER) of the model predictions f(X) and targets Y. For probabilistic models, the score func- tion (STD) is the absolute value of the 7-score (with the prediction mean, standard deviation and target given as fx))), fo (X), respectively. Having obtained distribu- tion of nonconformity scores & of the calibration dataset (X,, Yi) quantile I corresponding to the desired cov- erage 1 - is estimated from its cumulative distribution function F3 (Papadopoulos, 2008):
g'= = F-1 [nn 1)(1 a]]
(6)
The norm 1 (u)) of the residual operator itself pro- vides measure of UQ for the neural PDE. However, is limited by the accuracy of the gradient estimation method and can become computationally expensive when explor- ing a vast solution space (Tolsma & Barton, 1998). By using the residual norm as a nonconformity score, we fur- ther calibrate the approximate physics residual error that an inexpensive and coarse differential operator obtains. CP using PRE provides statistically valid and guaranteed error
n
The quantile estimates the error bar associated with the de- sired coverage and is combined with the new prediction to obtain the prediction sets. The nonconformity score func- tions and their prediction sets for AER and STD are given in Table 1. Both AER and STD are data-intensive CP methods, requiring calibration data that converges to a beta distribu- tion (Angelopoulos & Bates, 2023). This data dependence
3

---

Calibrated
UQ
bars across the PDE's residual space, incorporating physical information into the calibration procedure and providing calibrated measure of the physical misalignment of the surrogate model.
These coverage bands expand across the entire simulation domain 92 [0, 7] (discretised as NN,,NN,NN, rather than individual cell within it. For coverage band CO, the joint-CP formulation ensures that 1 - a predictions/solu- tions lie within the bounds. For performing joint-CP, the non-conformity scores are modified to reflect the supremum of the score functions in Table They are modulated by the standard deviation a of the calibration scores (Diquigio- vanni et al., 2021) to obtain prediction bands with varying widths based on local behaviour (Diquigiovanni et al., 2022). The modifications of the score functions and prediction sets to perform CP are given by
PRE as a nonconformity score enables data-free conformal prediction. The estimated scores rely only on the neural PDE predictions and not on the target as in AER and STD. The only criterion is that the calibration and prediction do- mains arise from exchangeable PDE conditions. As shown in Table 1, PRE gives "prediction sets" independent of the prediction inputs (see Appendix B for formalism). Tradi- tional CP methods rely on calibration using labelled data to construct prediction intervals that contain the true values with a specified confidence level a. These methods guaran- tee that the true solution will lie within the estimated error bounds based on this empirical calibration. Contrastingly, our CP-PRE formulation takes a fundamentally different approach. Instead of requiring target data for calibration, we leverage the unique property of PDEs where the true solution in the residual space exists at zero. This eliminates the need for empirical calibration data. Our method focuses on ensuring that predictions themselves fall within coverage bounds Co, rather than guaranteeing that the true solution lies within these bounds. This allows us to validate predic- tion sets without access to ground truth data- a significant advantage over traditional CP approaches. We formalise this novel property in our theoretical framework presented in Appendix A.
Syont sup Xen, (E
(8)
a
joint
= qth (8))
(9)
where á and p (1)++1) are the formulations of the nonconfor- mity scores and the prediction at Xx,+1 used for marginal- CP as shown in Table 1. Joint-CP becomes particularly useful in identifying predictions that fail to fall within cov- erage, allowing us to accept or reject prediction based on a predetermined probability. Similar to that demonstrated by Casella et al. (2004), our framework can perform acceptance- rejection using a CP-based criterion. The acceptance prob- ability based on confidence level a. If a prediction is rejected, the PDE parameters that led to those predictions could be provided to the expensive physics- bbsed numerical PDE solver for further evaluation as indicated in Figure 1.
4.1. Marginal-CP
4.3. Differential Operator: Convolutional Kernels
Stencils as
The CP formulation was initially conceptualised for cali- brating univariate functions with single-point outputs (Vovk et al., 2005). It has recently been extended to spatio- temporal data, with multi-dimensional outputs with an immutable tensor structure (Gopakumar et al., 2024a). Within such spatio-temporal settings, CP has been imple- mented to provide marginal coverage, i.e. the calibration procedure provides independent error bars for each cell within the spatio-temporal domain. For an output tensor Y E R. x X,, N,,,N where Nz, Ny. N₁ represent the spatio- temporal discretisation of the domain, marginal-CP uses the non-conformity scores outlined in Table 1 across each cell of Y to obtain error bars, which will be compliant with Equa- tion (5) for each cell. Marginal-CP using PRE helps indicate spatio-temporal regions within individual predictions that lie outside the calibrated bounds of physics violation and require specific attention, treating those prediction regions with caution.
Calibrating neural PDFs using PRE nonconformity scores requires frequent evaluations of the composite differential operator D in Equation (1). For PDEs, this involves es- timating spatio-temporal gradients across the discretised domain, ranging from millions in simple cases to billions of gradient operations for complex physics. To address this computational challenge, we developed scalable gradient estimation method for evaluating PRE.
We use convolution operations with Finite Difference (FD) stencils as convolutional kernels for gradient estimation (Actor et al., 2020; Chen et al., 2024a;b). For instance, the 2D Laplacian operator V2, using a central difference scheme with discretisation h, can be approximated by
0
1
0
1
-4
1
(10)
h2
0
1 0
4.2. Joint-CP
and used as a kernel. This approach is justified by the math- ematical equivalence of FD approximations and discrete convolutions. Both represent matrix- multiplications of a block Toeplitz matrix with a field vector (Strang, 1986;
The joint-CP formulation constructs a calibration proce- dure that provides coverage bands for multivariate functions.
4

---

Calibrated Physics-Inforiie UQ
Fiorentino & Serra, 1991). The efficiency of this method stems from the optimised implementation of convolution op- erations in machine learning libraries like PyTorch (Paszke et al., 2019) and TensorFlow (Abadi et al., 2015).
PDEs. We compare various Bayesian methods along with multivariate inductive conformal prediction our method as shown in Table 2. We demonstrate that our method is capable of providing valid coverage guarantees for both in and -it-o---iitriuution testing without a significant increase in evaluation times (including time taken for sampling in Bayesian methods, data generation in data-driven CP and subsequent calibration). Methods that attain the required coverage are emboldened. Table 2 provides a qualitative comparison of our method against the other benchmarks and highlights that our method is data-free, requires no mod- ification or sampling and provides guaranteed coverage in a informed manner. We confine our comparison studies in Table 3, 4, and to the Wave, Navier-Stokes and Magnetohydrodynamie equations.
The FD approximation offers several advantages over Auto- matic Differentiation (AD) for our application. It is compati- ble with CP as a post-hoc measure, requires no architectural modifications, and is model-agnostic. Furthermore, FD im- plemented via convolutions is more memory memory-efficient than AD, which requires storing the entire computational graph. As we focus on the (mis)alignment of neural PDEs with Equation (1), we disregard boundary conditions, but demon- strate that they can be accounted for with convolutional padding in Appendix (Alguacil et al., 2021). Utilising finite difference schemes, the residuals are estimated up to the truncation errors of the Taylor approximation (MacK- innon & Johnson, 1991). Although discretisation plays role in the width of the error bars, CP-PRE still guarantees coverage, and this is further explored in Appendix D.I.
The second campaign (Section 5.4, 5.5) applies CP-PRE to fusion applications. We enhance tokamak plasma be- haviour surrogate models to identify erroneous dispersion regions (Section 5.4) and integrate CP-PRE with tokamak design surrogates to identify viable designs and areas need- ing additional simulations (Section 5.5). This campaign demonstrates the utility of CP-PRE in complex, practical applications. One-dimensional PDE experiments demon- strating CP-PRE are demonstrated in Appendix F and can be reproduced using code in the supplementary material.
5. Experiments
W2ve
Navior
WHD
5.75
0.79
6.75
1.50
0.50
5.1. Wave Equation
0.25
6.25
0.23
The two-dimensional wave equation is given by
6:00
6.33
6.26
6.00
110
0.32
0.25
6.90
5.76
2.50
0.38
0.25
6.91
5.76
1.00
8u 0t2
82u
Pu.
Ideal
Mang jeent
(11)
On
dy²
Figure 3. Validation plots demonstrating coverage guarantee detailed in Equation (5) obtained by performing CP using PRE across experiments. The average empirical coverage obtained experimentally is given the y-axis (ranging from 0 to 1, with I representing 100% coverage), while the theoretical coverage is represented the æ-axis. We obtain guaranteed coverage while using marginal-CP formulation and near-to-ideal coverage for the Joint-CP formulation.
We solve the field u given in Equation (11) within the do- main x, y E [-1, t [0, 1.0] with wave velocity c= 1.0 using spectral solver with periodic boundary conditions (Canuto et al., 2007). The initial conditions are parame- terised by the amplitude and position of a Gaussian field. A 2D FNO is trained on this data to predict 20-time steps autoregressively from a given initial state. Figure com- pares the model predictions against ground truth, showing the PRE and confidence bounds from marginal and joint-CP at 90% coverage. The PRE reveals noise artefacts in regions where the field should be zero, highlighting physical incon- sistencies despite apparently accurate predictions. Joint-CP bounds are necessarily larger than marginal-CP bounds as they span across the spatio-temporal domain as opposed to being cell-wise. As demonstrated in Figure 3, both CP approaches achieve the expected coverage guarantees. Ta- ble 3 highlights the superior performance of CP-PRE over other methods of UQ where it guarantees coverage when evaluated for in and out-of--dstrrbutie without requiring any additional data, drastically reducing the evaluation time. Marginal-C3 -P shows linear coverage due to cell-wise averag-
CP-PRE experiments comprise two campaigns. First, we benchmark CP-PRE within standard neural PDEs (Sec- tion 5.1 to 5.3). The calibration process (Figure 2) involves: (a) sampling inputs to the model to generate predictions, (b) calculating PRE(s) scores, and (c) calibrating physical error using marginal and joint-CP formulations. Validation uses the same PDE condition bounds as calibration. Within this campaign, we compare our method (CP-PRE) with other methods of providing uncertainty estimation for neural-
'The Basic Linear Algebra Subroutines (BLAS) in these li- braries leverage vectorisation and efficient memory access for significant performance gains. Our experiments demonstrate 1000x speedup using frndk.nnfunnttoonl.. conv3d versus an equiva- lent numpy finite difference implementation on standard CPU.
5

---

Calibrated Physics-Inforin UQ
Method
Data-Free
Modification-Fee
Sampling-
Guaranteed Coverage
Dropout (Gal & Ghahramani, 2016)
X
X
X
Deep Ensemble (Lalshminarayanan et al., 2017)
X
X
BNN (MacKay, 1992)
X
X
XXX
SWA-G (Maddox et al., 2019)
CP-AER (Gopakumar al., 2024a)
X
CP-PRE (Ours)
Table Comparing features across various UQ measures. Our method is data-free, does not require any modifications or sampling, and
helps obtain guaranteed coverage bounds in innfrre- manner.
Prediction (u)
PRE (D(u))
Marginal CP (+q) 3
Joint CP (+q mod) x102
4.0
1.0
1.5
3.5
1.5
0.5
3.0
1.0
0.0
1.0
2.5
0055
2.0
0.5
1.5
-1.
Figure 4. Wave: (From left right) neural PDE (FNO) prediction the last time instance, physics residual error of the prediction, Upper error bars obtained by performing marginal-CP and joint-CP respectively (90% coverage). For brevity, we have only shown the upper bars of the symmetric prediction sets. mod represents the modulation function in Table 1. The physical inconsistencies within the residual space of the prediction are calibrated and bounded using CP-PRE.
ing. At the same time, joint-CP exhibits coverage variations that depend on the calibration dataset (see Table for sta- bility analysis across multiple calibration sets). Additional experimental details are provided Appendix J.
each equation. Table 4 shows that Bayesian methods fail to provide coverage for FNO modelling, while multivariate CP (CP-AER) provides coverage but requires substantially higher evaluation time arising from the generation of sim- ulation data for calibration. Figures and 21 show PRE bounds from both equations. Having two PDE residuals en- ables easier rejection of predictions violating both bounds. Physics details, parameterisation, and training (12) pendix K. Within the scope of this paper, we limit ourselves to measuring the deviation of the model with the PDE resid- ual. The CP-PRE formulation can be extended to obtain bounds for both the initial and boundary conditions, further explored within Appendix E.
5.2. Navier-Stokes Equation
Consider the two-dimensional Navier-Stokes equations
in Ap-
v = 0,
(Continuity equation)
(V VVV -
(13)
(Momentum equation)
5.3. Magnetolydrodynamies
where interested in modelling the evolution of the velocity vector (v v)) and pressure (P) field of incompressible fluid with kinematic viscosity (v). For data generation, Equations (12) and (13) are solved on a domain a E 1), E [0, 1), t E [0, 0.5) using a spectral-based solver (Canuto et al., 2007). A 2D multi-variable FNO (Gopakumar al., 2024b) is trained to model the evolution of velocity and pressure autoregressively up until the 20th time instance.
Consider the magnetohydrodynamic (MHD) equations
ap
(pV) = 0,
(Continuity equation)
(14)
Ot
1
Øt VV
= B x
B
VP.
(15)
140
(Momentum equation)
d
P
0,
(Energy equation)
(16)
Unlike the previous example, the Navier-Stokes case has two equations with corresponding PRE estimates: conti- nuity (Equation (12)) and momentum (Equation (13)) for mass and momentum conservation. Our CP-PRE method calibrates model deviation from physical ground truth for
d
p
ØB
x (V x B (Induction equation)
(17)
Øt
V. B 0,
(Gauß law for magnetism)
(18)
6

---

Calibrated Physics- UQ
in-distribution
oul-of-distribution
Time
UQ
MSE
Coverage (95%)
MSE
Coverage (95%)
Train (hr)
Eval (8)
Deterministic
1.77e-05 + 3.69e-07
2.46e-03 + 2.00e-05
0:38
22
MC Dropout
1.44e-04 4: 3.26e-06
97.31 + 0.03
2.12e-03 2.60e-05
89.83 0.07
0:52
120
Deep Ensemble
8.76c-06 2.43c-07
98.02 = 0.04
2.42c-03 = 1.58c-05
83.44 + 0.12
3:10
112
BNN
1.92e-04 1: 1.92e-06
97.10 1 0.09
2.67e-03 1.26e-05
91.76 0.10
0:53
118
SWA-G
1.4le-05 +: 1.74e-06
94.55 = 3.25
2.55e-03 2.82e-05
$1.90
0:47
113
CP-AER
1.76e-05 + 4.40c-07
95.70 + 0.21
2.46c-03 + 1.41c-05
95.59 + 0.14
0:38
2022
CP-PRE (Ours)
1.78e-05 + 4.61e-07
95.52 = 0.21
2.46e-03 + 1.25e-05
95.39 + 0.12
0:38
32
Table 3. Wave Equation - CP-PRE guarantees coverage across distributions while providing the quickest evaluation time.
in-distribution
out-of-distribution
Time
UQ
MSE
Coverage (93%)
MSE
Coverage (95%)
Train (hr)
Eval (s)
Deterministic
1.05e004 1 6.91e-06
3.67e-03 5.30e-05
3:22
25
MC Dropout
5.96c-04 + 2.30c-05
82.21 = 0.22
4.30c-03 + 8.05c-05
44.05 + 0.26
3:34
153
Deep Ensemble
1.22e-04 1: 3.95e-06
91.31 = 0.08
3.67e-03 3.52e-05
30.74 0.19
16:22
147
BNN
6.90e-03 + 1.31e-04
89.91 + 0.20
6.95c-03 + 1.31c-04
85.19 + 0.23
3:39
152
SWA-G
1.96e-04 + 1.15e-05
84.22 = 2.37
3.63e-03 + 1.37e-04
31.00 + 2.85
3:28
146
CP-AER
1.05e-04 + 6.58e-06
95.56 + 0.40
3.66e-033
2.81e-05
95.54 + 0.15
3:22
20026
CP-PRE (Ours)
1.07c-04 + 5.18c-06
95.44 = 0.22
3,70c-03 = 4.23e-05
95.57 = 0.14
3:22
134
Table Navier-Stokes Equations
CP-PRE guarantees coverage across distributions while providing the quickest evaluation time.
where the density (p), velocity vector (v = [u, uI) and the pressure of plasma is modelled under a magnetic field (B [B,, By]) across spatio-temporal domain z; y E [0, 1]2, e [0. pm is the magnetic permeability of free space. Equations (14) to (18) represent the ideal MHD equations as a combination of the Navier-Stokes equations for fluid flow with Maxwell's 's equations of electromagnetism (Alfvén, 1942; Gruber & Rappaz, 1985; Mocz et al., 2014). The equations assume perfect conductivity (no magnetic diffusivity) and no viscosity. We focus our experiment on the modelling of the Orszag-Tar Tang vortex of turbulent plasma (Orszag& & Tang, 1979) with the data being generated using a finite volume method (Eymard et al., 2000). A 2D FNO trained to model the evolution of all 6 variables over a dataset generated by parameterised initial conditions.
are significantly tighter, while joint bounds are wider but useful for identifying predictions that violate conservation equations. Additional CP plots using other residuals (Fig- ures 27 and 28) and physics/surrogate model details are in Appendix I..
5.4. Plasma Modelling within Tokamak
In (Gopakumar al., 2024b), the authors model the evo- lution of plasma blobs within fusion reactor (known as tokamak) using an FNO. They explore the case of electro- static modelling of reduced magnetohyydrodynamics with data obtained from the JOREK code (Hoelzl et al.,, 2021). In the absence of magnetic pressure to confine it, the plasma, driven by kinetic pressure, moves radially outward and col- lides with the wall of the reactor. The plasma is charac- terised by density p, electric potential & and Temperature T, and the FNO models their spatio-temporal evolution autore- gressively. Borrowing that pre-trained model and utilising the reduced-MHD equations within the toroidal domain, we demonstrate obtaining calibrated error bars using CP-PRE scale. The FNO demonstrated in (Gopakumar et al., 2024b) model the plasma six orders of magnitude faster than traditional numerical solvers, and by providing calibrated eтror bars over the predictions, wider range of plasma configurations can be validated.
Equations (14) (18) provide with five measures of estimating the PRE of the MHD surrogate model. Each diffferent set of variables as- sociated with the system, allowing us to infer errors con- tributed to each variable accordingly. Table demonstrates that the Bayesian methods fail miserably providing valid error bars over the MHD case, The CP-AER provides guar- anteed coverage but is heavily dependent on simulation data, adding the computational expense. CP-PRE using induction (Equation (17)) and energy (Equation (16)) are shown for 90% coverage (a 0.1), sliced at y 0.5m.
PRE estimate depends
We focus on the temperature equation within reduced-MHD (equation 3 within (Gopakumar et al., 2024b)) as it com- prises all the variables associated with the plasma. As shown
The plots show PRE along the x-axis with marginal and joint bounds at a specific simulation time, Marginal bounds

---

Calibrated Physics- UQ
in-distribution
oul-of-distribution
Time
UQ
MSE
Coverage (95%)
MSE
Coverage (95%)
Train (hr)
Eval (8)
Deterministic
2.20e-03 + 5.20e-03
4.71e-02 + 1.06e-03
5:00
40
MC Dropout
3.29e-02 +: 5.86e-04
41.13 + 0.19
2.09e-01 1.38e-03
16.91 0.06
5:30
240
Deep Ensemble
3.59c-03 3.51c-04
78.15 = 0.16
3,41c-01 3.15c-02
39.63 + 0.31
26:25
235
BNN
4.20e-03 1 4.08e-05
90.24 0.10
4.63e-02 8.98e-04
62.37 1 0.46
5:40
240
SWA-G
2.61e-03 i 9.68e-05
48.50 3.81
4.53e-02 6.64e-04
14.22 1.35
5:22
236
CP-AER
2.20e-03 + 4.38e005
95.61 + 0.26
4.69c-02 + 18c--004
95.60 + 0.27
5:00
40042
CP-PRE (Ours)
2.20e-03 + 4.96e-03
95.54 = 0.18
4.71e-02 1.06e-03
95.67 + 0.22
5:00
482
Table 5. Magnetohyydodynami Equations - CP-PRE guarantees coverage across distributions with only marginal increase in evaluation time, arising from residual evaluation over a larger family PDEs (see Equations (14) to (18)).
Absolute Error: (T)
PRE: D(p, 4.
where v/ represents the poloidal magnetic flux, the kinetic pressure, F rB the toroidal magnetic field, and Ho the magnetic permeability given in 2D toroidal coordinate system characterised by and While traditional numeri- cal solvers like EFIT++ and FreeGSNKE (Iaa al., 1985; Amorisco al., 2024) are used for equilibrium reconstruc- tion, their computational cost has motivated neural network alternatives (Joung et al., 2023; Jang et al., 2024). However, these surrogate models lack UQ capabilities, making their deployment within the control room challenging.
1
0
2
(a) Abs. Error
(b) PRE
Figure 5. Reduced MHD: CP-PRE using the Temperature equa- tion (Eqn. in (Gopakumar et al., 2024b)) of Adduced-MHll to bound the plasma surrogate models. The PRE captures the model error relatively well, allowing us to provide lower and upper error bars corresponding required coverage.
PRE: Das(V)
Lower
Upper
0.000
0.042
0.035
0.030
0.025
0.020
in Figure 5, our method can capture the model error across range of predictions and devise error bars that provide guaranteed coverage without additional data. In figure 5(a), we demonstrate the absolute егrоr in the model prediction of the temperature evolution, correlating that with the PRE over the temperature equations in figure 5(b). By obtaining bounds on the PRE, we can determine the efficacy of the surrogate model in evaluating plasma evolution and iden- tify conditions under which fails, running the MHD code JOREK under those failed conditions.
0.036
2.5 5.0
0.042
Figure 6. Grad-Shafranov: The PRE for a specific poloidal field coil configuration is indicated the left, and the lower and upper bars for 50% are displayed adjacent to it. Aside from guaranteeing coverage, the CP-PRE framework allows discard physically inconsistent equilibria predicted by the surrogate model.
5.5. Magnetic Equilibrium in Tokamak
We implement an auto-encoder that maps poloidal magnetic flux across the poloidal cross-section for given tokamak architectures, conditioned on poloidal field coil locations under constant plasma current. While this accelerates simu- lation by orders of magnitude, lacks physical guarantees. By incorporating Equation (19) within the CP-PRE frame- work, we identify physically stable equilibria and obtain statistically valid error bounds. Figure 6 shows the PRE over a surrogate model prediction with lower and upper er- ror bars for 50% coverage. Further details about the problem setting and the model can be found in Appendix N.
Tokamaks confine plasma within toroidal vessel using magnetic fields to achieve nuclear fusion. The plasma at high temperatures is contained by magnetic fields that coun- terbalance its kinetic pressure. This equilibrium state, function of magnetic coil configurations and plasma param- eters, is governed by the Grad-Shafranov (GS) equation (Somov, 2012)
86 1 1g
84
2 dp 1 dF
(19)
0r2 dr 0z2
du
2 dvp
8

---

Calibrated Physics-Informi UQ
6. Extensions of CP-PRE beyond PDEs
across multiple predictions (joint-CP). The latter enables rejection criterion for set of predictions, potentially serv- ing an active-learning pipeline for neural PDE solvers (Musekamp et al., 2025). CP-PRE provides guaranteed cov- erage irrespective of the model, chosen discretisation, or the PDE of interest; however, the width of the error bar indicates quantitative features in the model quality. well-trained model will exhibit tighter error bars as opposed to a poorer fit model; see Appendix I.
Our framework applies to any prediction case where the forward problem can be formulated as a residual with an equality constraint (Appendices A and B). The framework works in any scenario where the forward model can be expressed in the standard canonical form
Ar 0,
(20)
where is the model prediction, A is differential or al- gebraic operator governing the dynamics, and b is non- homogeneous term such as a function or a constant.
Our method's coverage bounds exist in the PDE residual space rather than the Euclidean space of phys- ical variables. Transforming to physical space involves chal- lenging set propagation through integral operations, which may require approximations (Teng et al., 2023) or expensive Monte Carlo sampling (Andrieu et al., 2003). The data-free approach lacks a grounding target for calibration, though we argue that large sample of model outputs provides sta- tistically significant overview of uncertainty. The sampling cost from the neural-PDE solver for calibration involves intensive gradient evaluations. PRE estimation using finite- difference stencils also introduces the errors associated with Taylor expansion. The current formulation is limited to regular grids with fixed spacing, though extensions un- structured grids via graph convolutions are possible (Eliasof & Treister, 2020). CP-PRE does not help differentiate be- tween aleatoric and epistemic uncertainty. It aligns with conformal prediction's 's characterisation of predictive un- certainty. From one perspective, this could be viewed as aleatoric uncertainty since we construct confidence intervals relative to a specific probability distribution (the distribu- tion from which calibration data, i.e. initial conditions/PDE coefficients, are sampled). Alternatively, it could be con- sidered epistemic uncertainty since we model the neural network's ss error through confidence intervals (typically used for unknown but fixed quantities). While we believe the latter interpretation is more appropriate, we acknowledge that the traditional aleatoric/epistemic dichotomy may not be directly applicable to our framework. This distinction is most valuable when both uncertainty types coexist and require separate treatment (Ferson & Ginzburg, 1996).
Limitations
This extends our applications beyond PDEs to ODEs and al- gebraic equations found in control problems (Jiang & Jiang, 2014), chemical reactions (Thöni et al., 2025), biological systems (Wang et al., 2019), and financial scenarios (Liu et al., 2019).
We deliberately focus on PDFs as they represent the most comprehensive and challenging ause-multi-dimensionee domains with complex spatio-temporal dependencies and unique computational challenges. Success with PDEs im- plicitly validates applicability to simpler systems.
7. Discussion
If "All models are wrong, but some are useful" (Box, 1976), through this work, we explore a novel framework for pro- viding a measure of usefulnesss of neural PDEs. We deploy a principled method of evaluating the accuracy of the solution, Le. its (calibrated) obedience to the known physics of the system under study. As opposed to other methods of UQ for neural PDEs, our method has the unique advantage of being physics-informed This allows us to study the physi- cal inconsistencies of the model predictions with coverage guarantees provided by conformal prediction. We conclude with a discussion of the strengths, limitations and potential improvements.
Strengths PRE estimates the violation of conservation laws in neural PDE predictions, guaranteed error bounds over the physics deviation. This post-hoc uncertainty quan- tification is model- and physics-agnostic scaling linearly with model complexity and quasi-linearly with PDE com- plexity due to the additive nature of differential operators. Our framework reformulates CP to be data-free, express- ing model inaccuracy solely through PRE, not requiring a labelled dataset. This approach reduces calibration costs and loosens exchangeability restrictions we can mod- ify the calibration and, hence, the prediction domain by simply reformulating the PRE accordingly. The PRE for- mulation (Section& 4, Appendix B) yields input-independee prediction sets, allowing for the identification of weak pre- diction regions within single simulations (marginal-CPP) and
8. Conclusion
We address the problem of the reliability of neural-PDD solvers by proposing CP-PRE, a novel conformal predic- tion framework. Our method provides guaranteed and physics-informin informed uncertainty estimates for each cell within prediction, identifying erroneous regions while discerning physically inconsistent predictions across the entire spatio- temporal domain. Our work enhances the reliability of neu- ral PDE solvers, potentially broadening their applicability in science and engineering domains where robust uncertainty quantification crucial.
9

---

Calibrated
UQ
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which feel must be specifically highlighted here.
10

---

Calibrated Physics-Informi UQ
References
0565. doi: 10. 10.1023//:10000081527716 URL https:
/ //dpi.org/10. 11tp://...0000111102222222111111
Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z.,
Citro, C., Corrado, G. S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowwi, R., Kaiser, L.., Kudlur, M., Lev- enberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M.. Shlens, J., Steiner. B., Sutskever, I. Talwar, K. Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., and Zheng. X. TensorFlow: Large- scale machine learming on heterogeneous systems, 2015. URL Utthttps: www.teensee. www.wllwwoooo org/ Software available from tensorflow.org.
Angelopoulos, A. N. and Bates, S. Conformal prediction:
A gentle introduction. Found. Trends Mach. Learn., 16 (4):494-591, mar 2023, ISSN 1935-8237. doi: 10. 10.1561/
2200000101. URL https: /twii/oiiio//// org/ 1561/ 2200000101.
Bartolucci, F., de Bezenac, E., Raonic, B., Molinaro, R., Mishra, S., and Alaifari, R. Representation equivalent neural operators: framework for alias-free operator learning. Thirty-seventh Conference on Neural In- formation Processing Systems, 2023. URL https: / //openzeview.co net forum?dd=
Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P., Cao, X., Khos- ravi, A., Acharya, U. R., Makarenkov, V., and Na- havandi, S. A review of uncertainty quantification in deep learning: Techniques, applications and chal- lenges. Information Fusion, 76:243-297 2021. ISSN 1566-2535. doi:
Bertone, G., Deisenroth, M. P., Kim, J. S., Liem, S., Ruiz de Austri, R., and Welling, M. Accelerating the bsm interpretation of lhe data with machine learning. Physics of the Dark Universe, 24: 100293, 2019.
Box, G. E. Science and statistics. Journal the American Statistical Association, 71 71(36))7791-79,, 1976. doi: 10. 1080/01621459.0 19w...176....999.
008. URLLAttppe//ww.seeeeneedccccccce com/ setiencearrticcepp////5555632331
Canuto, C.,, Hussaini, M., Quarteroni, A., and Zang, T. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. Scientific Com- putation. Springer Berlin Heidelberg, 2007. ISBN 9783540307280. URL https: /books..oss///o//o/oo google. co wo.uk/bae wuukkbooks7tao http/bbookkiid770gEw55....
Actor, J., Fuentes, D. T., and Riviere, B. Identification of kernels in a convolutional neural network: connections between the level set equation and deep learning for im- age segmentation. In Landman, B. A. and 1ggum, I. (eds.), Medical Imaging 2020: Image Processing. SPIE, March 2020.
Carey, N., Zanisi, L., Pamela, S., Gopakumar, V., Omotani, J., Buchanan, J., and Brandstetter, J. Data efficiency and long term prediction capabilities for neural operator surrogate models of core and edge plasma codes, 2024. URL https::///// /awxiv.vvvv. org/ abs / 08561........ 08561.
Alfvén, H. Existence of adecrooaagneii-lyydooyaaiii waves. Nature, 150(3805):405 406, Oct 1942. ISSN 1476-4687. doi: 10. 50.105000 doi. ori/o1g/10 1038 00088/1004550.
URL https://
Alguacil, A., Pinto, W. G., Bauerheim, M., Jacob, M. C., and Moreau, S. Effects of boundary conditions in fully dy-
Casella, G., Robert, C. P., and Wells, M. T. Gen- eralized accept-reject sampling schemes. Lec- ture Notes-Monograph Series, 45:342-344 2004. ISSN 07492170. URL http://www.ssorr.rrg.... stable/4356322
convolutional networks for learning spatio-tempont namics. In Dong, Y., Kourtellis, N., Hammer, B., and Lozano, J. A. (eds.), Machine Learning and Knowledge Discovery in Databases. Applied Data Science Track, pp- 102-117, Cham, 2021. Springer International Publishing. ISBN 978-3-030-86511--1
Chalapathi, N., Du, Y., and Krishnapriyan, A. S. Scal- ing physics-informed hard constraints with mixture- of-experts. In The Tweifth International Conference on Learning Representations, 2024. URL https:// openreview net forumm id=u3d2221322
Amorisco, N. C.. Agnello, A., Holt, G., Mars, M., Buchanan, J., and Pamela, S. Freegsnke: A python-based dy-
namic free-boundary toroidal plasma equilibrium solver. Physics of Plasmas, 3144)1042517, Apr 2024. ISSN 1070-664X. doi: 10.1063/5.01884677 URL https: //doi. wrw.ii.orr11..0 10w./006355... 0188467.
Chen, B., Heaney, C. E., Gomes, J. L., Matar, O. K., and Pain, C. C. Solving the discretised multi- phase flow equations with interface capturing on structured grids using machine learning libraries. Computer Methods in Applied Mechanics and En- gineering. 426:1 2024a. ISSN 0045-7825. doi:
Andrieu, C., de Freitas, N., Doucet, A., and Jordan, M. I. An introduction to meme for machine learning. Ma-
chine Learning, 50(1):5 Jan 2003. ISSN 1573-

---

Calibrated Physics- UQ
com/ sci snt.eee//arrrrce///////000553333
URL https://www.sieenee
multigrid pooling. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), Ad- vances in Neural Information Processing Systems, volume 33, Pp. 18016-18027. Curran Associates, Inc., 2020. URL http:t/////ooooeenng...eeeeee cc c//sapee_file/ppperr/////LL
Chen, B., Heaney, C. E., and Pain, C. C. Using ai libraries for incompressible computational fluid dynamics, 2024b. URL wwtpt::/tttttiieeeee
org/ abs/ 2402. 17913.
Crank, J. and Nicolson, P.
A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Mathematical Proceedings of the Cambridge Philosophical Sociery, 43(1):50-67, 1947. doi: 000077/00300041100000
d1 6a (aa277ccc66d00711229fffff0011113 Paper.. pdf.
Paper.
Eymard, R., Gallouët, T., and Herbin, R. Finite volume methods. In Lions, J. I., and Ciarlet, P. (eds.), Solution of Equation in Rn (Part 3), Techniques of Scientific Comput- ing (Part 3), volume 7 of Handbook of Numerical Analy- sis, pp. 713-1020. Elsevier, 2000. ISBN 9780444503503. doi:
Danabasoglu, G., Lamarque, J.-F., Bacmeister, J., Bailey, D. A., Duviee,r, A. K., Edwards, J., Emmons, L. K., Fasullo, J., Garcia, R., Gettelman, A., Hannay, C., Hol- land, M. M., Large, W. G., Lauritzen, P. H., Lawrence, D. M., Lenaerts, J. T. M., Lindsay, K., Lipscomb, W. H., Mills, M. Neale, R., Oleson, K. W., Otto-Bliesner, B., Phillips, A. S., Sacks, W., Tilmes, S., van Kampenhout, L., Vertenstein, M., Bertini, A., Dennis, J., Deser, C., Fischer, C., Fox- Kemper, B., Kay, E., Kinnison, D., Kushner, J., Larson, V. E., Long, M. C., Mickelson, S., Moore, J. K., Nienhouse, E., Polvani, L., Rasch, P. J., and Strand, W. G. The community earth system model version 2 (cesm2). Journal of Advances in Modeling Earth Systems, 1202):22019MS000916, 2020. doi: http:://ooi.oog//100099991900009
Ferson, S. and Ginzburg. L.. R. Different methods are needed to propagate ignorance and variabil- ity. Reliability Engineering & System Safery, 54(2):133-144, 1996. ISSN 0951-8320. doi:
URL http:://www.ccieneed.ree. com/ ccckenee/srttcdeee///////e122222
Treatment of Aleatory and Epistemic Uncertainty.
Fiorentino, G. and Serra, S. Multigrid methods for toeplitz
matrices. CALCOLO, 28(3):283-305, Sep 1991. ISSN 1126-5434. doi: 10.1007/BF025758866 URL https: //dpi. 10 1007 B000000///2255888888
URL
https:://gguaaa wwwp://wwwoo..nnnnnnneeneneaamee wiley.
on
homm//oo/abs///..oe 9000000090000099009911 c2019MS001916 2019MS001916. 900100000999100000000009....
Gal, and Ghahramani, Dropout as a bayesian approxi- mation: Representing model uncertainty in deep learning. In International Conference Machine Learning, 2016.
de Hoop, M. V., Huang, Z., Qian, E., and Stuart, A. M. The cost-accuracy trade-off in operator learning with neural networks, 2022. URL //appsv//aaxvv.../ abs/2203.12 13181.
Geneva, N. and Zabaras, N. Modeling the dynam- ies of pde systems with physics-constrainee deep auto-regressive networks. Journal of Computa- tional Physics, 403:109056, 2020. ISSN 0021- 9991. doi: URL (tteencc/arrccle////////29999999
Diquigiovanni, J., Fontana, M., and Vantini, S. The impor- tance of being a band: Finite-sample exact distribution- free prediction sets for functional data, 2021. URL https: http///aavvvv/ee///e//eeeeeaaeva org/ abs/ 2ww2.112220 6744.
Diquigiovanni, J., Fontana, M., and Vantini, S. Conformal prediction bands for multivariate functional data. Journal of Multivariate Analysis, 189: 104879, 2022. ISSN 0047- 259X. doi: sttpp::/www.scienceet com/ URL sci sccience/artiicle////////0333333
Giles, D., Briant, J., Morcrette, C. J., and Guillas, S. Em- bedding machine-learn sub-grid variability improves climate model precipitation patterns. Communications Earth Environment, 5(1):712, Nov 2024. ISSN 2662- 4435. doi: 10. 10.10388/44227-004400888888 URL https: 01885-8-8-
Giudicelli, G., Lindsay, A., Harbour, L., Icenhour, C., Li, M., Hansel, J. E., German, P., Behne, P., Marin, O., Stogner, R. H., Miller, J. M., Schwen, D., Wang, Y., Munday, L., Schunert, S., Spencer, B. W., Yushu, D.,. Recuero, A., Prince, Z. M.. Nezdyur, M., Hu, T., Miao, Y., Jung, Y. S., Matthews, C., Novak, A., Langley, B., Truster, T., Nobre, N., Alger, B., Andrš, D., Kong, F,, Carlsen, R., Slaughter, A. E.,
Du, Y., Chalapathi, N., and Krishnapriyan, A. S. Neu- ral spectral methods: Self-supervised learning in the spectral domain. In The Twelfth International Confer- ence on Learning Representations, 2024. URL https:
Eliasof, M. and Treister, E. Diffgen: Graph convolutional networks via differential operators and algebraic
12

---

Calibrated
UQ
Peterson, J. W., Gaston, D., and Permann, C. 3.0 MOOSE: Enabling massively parallel multiphysies simulations. SoftwareX, 26:101690, 2024, ISSN 2352- URL.
Holod, I., Hu, D., Kim, S. K., Korving, S. Q.. Kos, L., Krebs, I., Kripner, L., Latu, G., Liu, E., Merkel, P., Meshcheriakov, D., Mitterauer, V., Mochalskyy, S., Morales, J. A., Nies, R., Nikulsin, N., Orain, F. Pratt, J., Ramasamy, R., Ramet, P., Reux, C.,, Särkimäki, K., Schwarz, N., Verma, P. S., Smith, S. F. Sommariva, C., Strumberger, E., van Vugt, D. C., Verbeek, M., West- erhof. E., Wieschollek, F., and Zielinski, The jorek non-linear extended mhd code and applications to large- scale instabilities and their control in magnetically con- fined fusion plasmas. Nuclear Fusion, 61 61(6))06500, 2021. doi: 10.1088//441-4456/a1199i //dx. doi www.ooggg/o...o
7110. doi: httpe://wwwsscceeeeeiireeeeeree c
Gopakumar, V. and Samaddar, D. Image mapping the tem- poral evolution of edge characteristics in tokamaks using neural networks. Machine Learning: Science and Tech- nology, ((1):015006, 2020. doi: 10.1008//232-21553/ ab5639. URL https: / /xx. doi wttp:oooi.oo...110000 1152-22533ab533
URL https: aba599..
Gopakumar, V., Pamela, S., and Samaddar, D. Loss landscape engineering via data regulation on PINNs. Machine Learning with Applica- tions, 12:100464, 2023a. ISSN 2666-8270.
Horwitz, J. A. K. Estimating the carbon footprint of compu-
tational fluid dynamics. Physics of Fluids, 36(4):045109 04 2024. ISSN 1070-6631. doi: 10. 1063/00030000 00093/5//00095.... URL htwps:tpps/ooom.....m doi. www.oo.g//... 063 / 5. 01 99350.
doi: http:///bbiorg//000016.....m com/ URL https: https://wwwwscccceceeeerettt /riance/Articcee/ii///6666683002
Iserles, A. A first course in the numerical analysis of differ- ential equations. Cambridge university press, 2009.
Jang, B., Kaptanoglu, A. A., Gaur, R., Pan, S., Lan- dreman, M.. and Dorland, W. Grad-shafranov equi- libria via data-free physics informed neural networks. Physics of Plasmas, 31(3):032500, Mar 2024. ISSN 1070-664X. doi: 10.1063/5.018866 URL https: ///oi.oor///// 00tg//...163/5... 0188634.
Gopakumar, V., Pamela, S., Zanisi, L.. Li, Z., Anandkumar, A., and Team, M. Fourier neural operator for plasma mod- elling, 2023b. URL https::///......... /araiv/////ivvvo//aa org/ abs/ 2602.0
Gopakumar, V., Gray, A., Oskarsson, J., Zanisi, L.., Pamela, S., Giles, D., Kusner, M., and Deisenroth, M. P. Uncer- tainty quantification of surrogate models using confor- mal prediction, 2024a. URL https: http://axxvvvo/// org/ abs/2408 09881.
Jiang, C. M., Esmaeilzadeh, S., Azizzadenesheli, K., Kashinath, K., Mustafa, M., Tchelepi, H. A., Marcus, P., Prabhat, and Anandkumar, A. Meshfreeflownet: A physics deep continuous space- time super- resolution framework. arXiv:2005.01 2020.
Gopakumar, V., Pamela, S., Zanisi, L.,, Li, Z., Gray, A., Brennand, D., Bhatia, N., Stathopoulos, G., Kusner, M., Deisenroth, M. P., Anandkumar, A., the JOREK Team, and Team, M. Plasma surrogate modelling using fourier neural operators. Nuclear Fusion, 64(5):056025, 4 2024b. doi: 10. 00.1088//74411442244d1111111 doi.orrg11 11/ 1088100880000 114111111111 431111433//ddd31
Jiang, Y. and Jiang, Z. P. Robust adaptive dynamic program- ming and feedback stabilization of nonlinear systems. IEEE Transactions on Neural Networks and Learning Systems, 25(5):882-893, 2014. doi: 10.1109/TNNLS. 2013.2294968.
URL https:////// //dx.
Jordà, M., Valero-Lara, P., and Peña, A. J. Performance evaluation of cudnn convolution algorithms nvidia volta gpus. IEEE Access, 7:70461-70473, 2019. doi: 10. 10//A/CESSS
Gray, A., Gopakumar, V., Rousseau, S., and Destercke, S. Guaranteed confidence-band enclosures for pde sur- rogates, 2025. URL https: /ttts::/arxiv... org/ /bs/ 2501.18426.
Joung, S., Ghim, Y.-C., Kim, J., Kwak, S., Kwon, D., Sung, C., Kim, D., Kim, H.-S., Bak, J. G., and Yoon, S. W. Gs-deepnet: mastering tokamak plasma equilibria with deep neural networks and the grad-shafranov equation. Scientific Reports, 13(1):15799 Sep 2023. ISSN 2045-
Gruber, R. and Rappaz, J. The Ideal MHD Model, pp. 34-11. Springer Berlin Heidelberg, Berlin, Heidel- berg, 1985. ISBN 978-3-642-86708-8,
doi: 10.1007/ URL httttts::/dii.oq////// 0111//98883-4-2--
2322. doi: 100000//44119-002222999999999999 URL https: htt.:/1..1101//111 1tw.1110///41--0000 //2.:..rr///0.000///4319--0----0
Hoelzl, M., Huijsmans, G. T. A., Pamela, S. J. P., Bécoulet, M., Nardon, E., Artola, F., Nkonga, B., Atanasiu, C. V., Bandaru, V., Bhole, A., Bonfiglio, D., Cathey, A., Czarny, O., Dvomova, A., Fehér, T., Fil, A., Franck, E., Fu- tatani, S., Gruca, M., Guillard, H., Haverkort, W.,
Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L. Physics-informed ma- chine learning. Nature Reviews Physics, 3(6):222- 440, 2021. ISSN 2522-5820. doi: 10.1038/
13

---

Calibrated Physics-Informi UQ
$42254-021-0033145. https: httpp/::///oi.oo/o/.//// 1038/$42254-021-03333
Liu, S., Oosterlee, C. W., and Bohte, M. Pricing op- tions and computing implied volatilities using neural networks. Risks, 7(1,, 2019, ISSN 2227-9091. doi: 10.3390/riss:7010016 URL attps://wwt https://www.mdi....
Kato, Y., Tax, D. M., and Loog, M. A review of noncon- formity measures for conformal prediction in regression. In Papadopoulos, H., Nguyen, K. A., Boström, H., and Carlsson, L. (eds.), Proceedings of the Twelfth Sympo- sium on Conformal and Probabilistic Prediction with Applications, volume 204 of Proceedings of Machine Learning Research, pp. 369-383, PMLR, 13-15 Sep 2023. https::: / /ptp/ceeeddngs...///////////..se mlr press// 2204/kato200 html.
Ma. Z., Azizzadeneshelii K., and Anandkumar, A. Cali- brated uncertainty quantification for operator learning via conformal prediction. arXiv preprint arXiv:2402.011 40200166, 2024.
MacKay, D. J. C. A practical bayesian framework for back- propagation networks. Neural Computation, 4(3):448- 472, 05 1992. ISSN 0899-7667. doi: 10.1162/neco.
Kingma, D. P. and Ba, J. Adam: A method for stochas- tic optimization. In Bengio, Y. and LeCun, Y. (eds.), 3rd International Conference on Learning Represen- tations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http: // //arxiv. org/ abs//411111 6980.
1992.4.3.448. URL https: /tctpps://doo.. /ttp:/:oooooooo/.// 1262/ neco 19900.192..0 1.3.1.999.......... 4.3.448..
MacKinnon, R. J. and Johnson, R. W. Differential- equation-based representation of truncation errors for accurate numerical simulation. International Journal for Numerical Methods in Fluids, 13(6):739-757) 1991. doi:
Kochkov, D., Yuval, J., Langmore, I., Norgaard, P., Smith, J., Mooers, G., Klöwer, M., Lottes, J., Rasp, S., Düben, P., Hatfield, S., Battaglia, P., Sanchez-Gonzalez, A., Will- son, M., Brenner, M. P., and Hoyer, S. Neural gen- eral circulation models for weather and climate. Na- ture, 632(8007) 4687. doi: 10. 10.10038/4415586-02--777444444 URL https: 11w.1118////30 wttp:////....1//111.e...
16001 1650130606.
URL eey. cwm/coooooooooooo
Aug 2024. ISSN 1476-
Maddox, W. J., Izmailov, P., Garipov, T., Vetrov, D. P., and Wilson, A. G. A simple baseline for bayesian uncertainty in deep learning. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché- Buc, F., Fox, E., and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.,
Lakshminarayanan, B., Pritzel, A., and Blundell, C. Simple and scalable predictive uncertainty estimation using deep ensembles, 2017.
2019. URL https://rrccce ngs neerrippp. cc/paper cc/sappe__iiee/ppppr///////LL 1 11B921efa223c3226666027886662222 pdf.
Lao, L., St. John, H., Stambaugh, R., Kellman, A., and Pfeiffer, W. Reconstruction of current profile parame- ters and plasma shapes in tokamaks. Nuclear Fusion, 1611 nov 1985. doi: 10.1088/0029-5515/33) 11/007. URL https: //dxxddo org.io..g100...000. 0029 51155555151//////1
McCabe, M., Harrington, P., Subramanian, S., and Brown, J. Towards stability of autoregressive neural opera- tors. Transactions on Machine Learning Research, 2023, ISSN 2835-8856. URL https: /apps////peeeeiee..e net/ forum? d- FOUtKYOG.
Li, S., Jiang, H., Ren, Z.,, and Xu, C. Optimal tracking for divergent-type parabolic pde system in current profile control. Abstract and Applied Analysis, 2014:1-8, 06 2014. doi: 10. 00.1155/0144940966
Meyer, D., Hogan, R. J., Dueben, P. D., and Mason, S. I. Machine learning emulation of 3d cloud radiative effects. Journal of Advances in Model- ing Earth Systems, 14)))22021 2022. doi: http:///bo..rr//000022222222222 https://aapguubb..ooiieeee wile.. c2021MS002550 2021M2002550
Li, Z., Kovachki, N. B., Azizzadenesheli, K., liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Fourier neu- ral operator for parametric partial diffcrential equations. In International Conference on Learning Representations, forum? fww.uuuummmet//uu
URL
2021. URL https: / http///wwwwwe..eeeeoppeeoeoeeeet net id=c8P 9N0Vtmno.
Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., and Anandkumar, A. Physics- informed neural operator for learning partial differential equations. ACM/ / IMS Data Sci., 1(3), may 2024. doi: 10. 10.1145/3648566 URL https: /ttp://///o........./ooooo/o 1145/364855
Meyer, H. and Team, S. Plasma burn mind-miidd-ii the gap. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 382 (2280):00300466, 2024. doi: 10.10998/ssaa.2223..000 URL http::://ooaalooceerrpppaaaaeae org/dol/abs/10. 2023..0....2222200...00
14

---

Calibrated Physics-Informen UQ
Moez, P., Vogelsberger, M., and Hernquist, L. A con- strained transport scheme for MHD unstructured static and moving meshes. Monthly Notices of the Royal Astronomical Sociery, 442(1):43-55, 06 2014. ISSN 0035-8711. doi: (0.10///mars///uu8865. URL https: //doi. 10. 1093
Raissi, M., Perdikaris, P., and Karniadakis, G. Physics- informed neural networks: A deep learning frame- work for solving forward and inverse problems involv- ing nonlinear partial differential equations. Journal of Computational Physics, 378:686-777, 2019. ISSN 0021-9991. doi: https///oiiiorg//00000/00/222220
045. URL htt.htsss/www.ccciceeeeeeee /tbeee///atiidee//////0299999999
com/
Musekamp, D., Kalimuthu, M., Holzmüüler, D., Takamoto, M.. and Niepert, M. Active learning for neural PDE solvers. In The Thirteenth International Conference Learning Representations, 2025. URL https: openreview net/ forumoirrriooo 11= Zm@aumRg.
Saad, Y. and Schultz, M. H. Gmres: A generalized mini- mal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Com- puting, 833) 1986. doi: 10. 1137/0907058. URL https: /tt////dii org/ 111100..777770055
Orszag, S. A. and Tang, C.-M. Small-scale structure of two- dimensional magnetohydrodynamic turbulence. Journal of Fluid Mechanics, ()9 1979. doi: 10.1017/ $00221 12079001XXX
Shafer, G. and Vovk, V. A tutorial on conformal prediction. Journal of Machine Learning Research, 9(3), 2008
Papadopoulos, H. Inductive conformal prediction: Theory and application neural networks. In Fritzsche, P. (ed.), Tools in Artificial Intelligence, chapter 18. IntechOpen, Rijeka, 2008. doi: 10.5772/6078.
Somov, B. V. Plasma Equilibrium in Magnetic Field, pp. 403-427. Springer New York, New York, NY, 2012. ISBN 178-1-4614-4228-- 978-1-461444283-7.199 URL /ttts:://oi.... 10 10.1007/97-1144144-111
doi: 10.1007/
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., Vitioo, Z., Rai- son, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. Pytorch: An imperative style, high-performance deep learning library, 2019. URL https: /ttp///axivv....///iv/e../.....a org/ wtwp///g////ooo.eee..e/eeeeeebee /bbs///11..110
Stankeviciute, K., M. Alaa, A., and van der Schaar, M. Con- formal time-series forecasting. In Ranzato, M., Beygelz- imer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), Advances in Neural Information Processing Sys- tems, volume 34, pp. 6216-6228, Curran Associates, Inc., 2021. URL https: /prope/d///cccediggaa neurips.... cc cccpppppe__ileeeepppee/////////2 31 2fffbaa2a723..eedaae95553333 pdf.
Pfaff, T., Fortunato, M. Sanchez-Gonzalez, A., and Battaglia, P. Learning mesh-based simulation with graph networks. In International Conference on Learning Rep- resentations, 2021. URL https:: //ppnreeveeeeseeeeeeee net net//Eorum2id-roNqyY0
Paper.
Strang. G. A proposal for toeplitz. matrix calculations.
Studies in Applied Mathematics, 74(2):171-176, 1986. doi:
URL
Pinder, G. F. Numerical methods for solving partial dif- ferential equations: : comprehensive introduction for scientists and engineers. John Wiley and Sons, Inc. Wiley, Hoboken, NJ, 2018.
wwwp:iiwwwa..e.....e.o...mooeeee com/doi/ abs/10. 0ttp:/00.10022/amm11186444444
Sun, S. Conformal methods for quantifying uncertainty in spatiotemporal data: A survey, 2022.
Poels, Y., Derks, G., Westerhof, E., Minartz, K., Wiesen, S., and Menkovski, V. Fast dynamic ld simulation of divertor plasmas with neural pde surrogates. Nuclear Fu- sion, 63(12):126012 sep 2023, doi: 10. 10.8////17414432/ org/10. 1088/
Teng, J., Wen, C., Zhang, D., Bengio, Y., Gao, Y., and Yuan, Y. Predictive inference with feature conformal prediction. In The Fleventh International Conference on Learning Representations, 2023. URL https: btttp:://openeeviert net / Sarrmmmdd=0ummmm YmFTu.
acf70d. URL https: 1741 //326/acf77
Thöni, A. C. M., Robinson, W. E., Bachrach, Y., Huck, W. T. S., and Kachman, T. Modelling chemical reaction net- works using neural ordinary differential equations, 2025. URL https: arwiv org/ abs 221///202.1133
Psaros, A. F. Meng, X., Zou, Z., Guo, L., and Karniadakis, E. Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons. Journal of Computational Physics, 477:111902, 2023. ISSN 0021-
9991. doi: https://ooiaagg//00.00.jg.gppp.2 URL bttp:://www.scienccddeeeett coclanee//rticee////////22292220
Tolsma, J. E. and Barton, P. I. On computational differentiation. Computers and Chemical Engi- neering, 22(4):475 1998. ISSN 0098-1354.
15

---

Calibrated Physics-Inforiit UQ
doi: http://dh..a///...oooss///851110
URL stt.:::/www..ccieneeeiiiie ccceneesarrrcee//////00111111111
sci
Vovk, V., Gammerman, A., and Shafer, G. Algorithmic Learning Random World. Springer, 2005.
Wang, S., Fan, K., Luo, N., et al. Massive com- putational acceleration by using neural networks to emulate mechanism-based biological models. Na- ture Communicationss 10:4354, 2019. doi: 10.1038/ $41467-019-12342-y, URL https: /ttp:/oo...oo../////////o
Xu, C. and Xie, Y. Conformal prediction interval for dy- namic time-series. In Meila, M. and Zhang, T. (eds.), Proceedings of the 38th International Conference on Ma- chine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 11559-11550 1569.. PMLR, 18-24 Jul 2021. URL https: ststtps://proceedinng mlr / v139 //wplhhhh.. html.
Xu, C., Xie, Y., Vazquez, D. A. Z.,, Yao, R., and Qiu, F. Spatio-temporal wildfire prediction using multi-modal data. IEEE Journal on Selected Areas in Information Theory, 4:302-313, 2023. doi: 10.1109//SATT.2223. 3276054.
Zhu, Y., Zabaras, N., Koutsourelakis, P.-.., and Perdikaris. P. coostreec deep learning for high-dimensional surrogate modeling and uncer- tainty quantification without labeled data. Journal of Computational Physics, 394:56-81, 2019. ISSN 0021-9991. doi: http:///oo..ay//o....6.jg.g.2pp9
024. URL U.t.etppss/.ww..ccceeeeeeeeeeeee ociianec/arrticlle///////2299999
Zou, Z., Meng, X., Psaros, A. F., and Karniadakis, G. E.
Neuraluq: comprehensive library for uncertainty quan- tification in neural differential equations and operators.
SIAM Review, 66(1):161-190, 2024. doi: 10.1137/
22M1518189 199 URL https:::::/ /ttp///oooo.....//////// 1137/
22M1 518189.
16

---

Calibrated Physics-Informet UQ
A. Theorem: Data-Free CP
B. PRE: Score Function and Prediction Sets
Preliminaries: Let D: R + Rm be a physics residual operator mapping function to its PDE residual value, where: {x is the calibration set, f is the model, qe is estimated the [(n 1)(1 a)l/n -quantile of {DD(f(X)))
For a general nonconformity score S, the prediction set for new input Xn+1 is typically defined as:
("(Xn+1) = (y ((x +)
where q is the (1 - a -quantile of the nonconformity scores on the calibration set.
Theorem A.1. If the residuals {D(f(X)))EL are ex- changeable random variables, then for any significance level (0, 1) and any new input Xu+1 have the fol- lowing coverage guarantee:
For AER and STD, the nonconformity scores depend on both the input X and the output (target) Y:
SAFR(X,Y) 1f(x) - YI,
(21)
P(I()((n+1)) ()) 0.)21-a; Co = -1- dn.
SsTD(X, (fu(X)- -VI
(22)
fa(X)
Proof. Let R ID(f(X.)) for I..... + We have, by assumption, (Ri,..., Rys Ra41) is exchange- able sequence. Define the rank of Ru+1 w.r.t. all other residuals:
The resulting prediction sets are:
CAER(Xn+1) = [f(Xn+1) - 8¹ ((xn11) + +"), (23)
T ((n+11) (i=1, +1:+1. Ratill
CSrD(Xnt1) = (fu(Xn11) 8 fo (Xn+1))
f (X fo (Xn+1))
(24)
By exchangeability, the rank ot (R)+1) is uniformly dis- tributed over (1...., 13. Therefore,
These prediction sets clearly depend on the input X,+1-
((T(Bn+1) < [())))(-a))) [[n + 1)(1-a)) aa]
For PRE, the nonconformity score depends only the model output and not the target:
2 1-a.
By construction of g we have that,
SprE(f(X) ID(((())) 0/,
{T(Rn+1) S [(++ 1)(1-)00)))))) a)]3 C (Ba+1 S }.
where D is the PDE residual operator. The key difference is that the true output Y for PRE, irrespective of the PDE is always O and does not depend on the input X. PRE is measure of how well the model output satisfies the physics rather than how it fits certain data. Hence, we can formulate a nonconformity score that is data-free and eventually leads to input-independent prediction sets as given below.
Putting this together,
P(ID(f(X s²*)= P(R 199)
which completes the proof.
For PRE, reframe the prediction set definition:
CT (f(x): (P(((X))) 4")
This set is not defined in terms the true values but in terms of the allowable model outputs f(X) that satisfy the PDE residual constraint. Thus, the prediction set can be expressed as:
CpRE = 1-
This formulation is independent of the input X, as it only depends on the quantile gDt derived from the calibration set as given in Equation (6).
To validate predictions using PRE:
For a new input Xn+1, compute f(Xn+1)-
Calculate the residual: r r= r=|D(((X,,+)))
17

---

Calibrated Physics-Informer UQ
3. Check if [-1", [ for given
(a) For marginal conformal prediction: Apply steps
1-5 cell-wise across the spatial-temporal domain
If the condition in step 3 is satisfied, the error bounds dic- tated by 1-go, (0) is considered valid according to the CP framework, regardless of the specific input X,+1-
to get localized error bounds
(b) For joint conformal prediction: Calculate the supremum the normalized residuals across the entire domain obtain global error bounds
C. Algorithmic Procedure
(c) The width of the error bounds indicates the model's ' physical consistency tighter bounds sug- gest better alignment with the physics
1. Set up the Neural PDE Solver
(a) Define the PDE system of interest with its govern- ing equations in numerical solver
7. Decision Framework (Optional)
(b) Train neural network (e.g., Fourier Neural Oper-
(a) Accept predictions where the physics residual falls within the calibrated bounds
ator) approximate solutions to the PDE
(c) Ensure the model can make predictions on new initial conditions/ PDE coefficients
(b) Reject and flag predictions where the residual ex- ceeds the bounds, potentially routing these cases to traditional numerical solvers
2. Define the Physics Residual Error (PRE)
(c) This approach guarantees that, with probability least a, the physics residual of new predic- tions will fall within the computed bounds, ensur- ing that the model maintains physical consistency while providing faster solutions than traditional numerical methods.
(a) For a given PDE case with operator D, the physics residual defined as [Df(())), where f is the neural PDE prediction
(b) Create differential operator using finite differ- ence stencils implemented as convolutional ker- nels. This operator evaluates how well the pre- dicted solution satisfies the underlying physics
3. Generate Calibration Set
(a) Sample set of initial conditions X1, X2.... Xn from the domain of interest
(b) Run the neural PDE solver on these initial condi- tions to get predictions
(c) Compute the PRE nonconformity scores for
DD(((x))), (I,,XX)))...,II((((X,, each prediction
4. Calibration Process
(a) For a desired confidence level 1-a (e.g., 90% for
a = 0.1), compute the empirical quantile:
(b) g
[(n + 1)(1 /n quantile of D()(X))) ((f(a)))), (D(,,,,))(((X)))
(c) This quantile represents the threshold for the physics residual that will provide the desired cov- erage
5. Apply to New Predictions
(a) For new initial condition X,+1,, obtain the neu- ral PDE prediction f(X,+1)
(b) Compute the physics residual [D(f(Xn+1))
(c) The prediction set defined as C -- (e)
(d) If (X(((x+1)) | falls within Co the prediction is considered valid at the 1 confidence level
6. Interpretation of Results
18

---

Calibrated Physics-Informed UQ
D. ConvOperator: Convolutional Kernels for Gradient Estimation
alpha, bsta 1.0, 0.5 Example
coefficients
D ConvOperator() () #Additive Kernels D.kernll D Det.eernel - alpha xx
Within the code base for this paper, release utility function that constructs convolutional layers for gradient estimation based on your choice of order differentiation and Taylor approximation. This allows for the PRE score function to be easily expressed in a single line of code
.kernel beta D i_deettity.keree
The convolutional kernels are additive i.e. in order to esti- mate the residual in one convolutional operation, they could added together to form a composite kernel that charac- terises the entire PDE residual.
This section provides an overview of the code implementa- tion and algorithm for estimating the PRE using Convolution operations. We'll use arbitrary PDE example with tem- poral gradient ăr the process.
Once having up the kernels, PRE estimation is as simple passing the composite class instance D the predictions from the neural PDE surroga te (ensuring that the output is in the same order as the kernel outlined above).
8y and a Laplacian
8
to illustrate
Oy
y P_pred model (X)
Ou
8 2x Oy²
8u
D (y (y_prod)
+ Bu -0,
(25)
Øt
Only operating on the outputs, this method of PRE estima- tion is memory efficient, computationally cheap and with the ConvOperator evaluating the PDE residual can be done in a single line of code.
where is the field variable, is time, and y spatial coordinates, and and are constants. To estimate the PDE residual given by Equation (25), we need to estimate the associated spatio-temporal gradients.
First, we use the ConvOperator class from set up the convolu-
Utils /ossls/convQps.dd.yy tional layer with kernels taken from the appropriate finite difference stencils:
from ConvÖps 2d import ConvOperator
Define each operator within the PDE D_t = ConvOperator (domain='t', crder=1) #time-derivative
D_xxyy CenvOperator (domain=('x","
order=2) #Laplacian
D b_identity ConvOperator () #Identity
Operator
The ConvOperator class is used to set up gradient operation. It takes in variable (s) of differentiation
order of di ifeerntiation as arguments to design the ap- propriate forward difference stencil and then sets up convolutional layer with the stencil as the kernel. Under the hood, the class will take care of devising convolutional layer, and setup the kernel so that it acts on spatio-temporal tensor of dimensionality: (BS, Nt, Nx, Ny] which expands to batch size, temporal discretisation and the spatial discretisation in x and y.
and
3D
"The code and associated utility functions be found in: https: / hggttubb../tggggggg///ggg....h.. wwwc www.cc/.cciieeeeeeeeeeeeeeeee..e gitvicky/ccsttteeeecccccc PriikkyCC-PPRS
"Convolution functions can be set as cross-correlan cross-oorrlltions as it default in the PyTorch framework, they could be set up convolutions by obtaining the complex conjugate in the frequency space.(Jorda) et al., 2019)
19

---

Calibrated Physics-Informue UQ
D.I. Impact of Discretisation
As demonstrated in (Bartolucci al., 2023), the discretisa- tion of the imputs and hence model outputs plays impor- tant role in the accuracy of the neural-PDE solvers. Though the neural operators are constructed for discretisation- invariant behaviour due to the band-limited nature of the be-
Coarse [dt=0.01)
0.015
PRE Upper Bound
0.010
functions, they often exhibit dlbccrttisationccovvrrr haviour as opposed to discretisation- -nvaaiiineene.....e...... This is of particular importance in the temporal dimensions as these neural-PIEE models utilise discrete, autoregressive based time-stepping, baked into the model within its train- ing regime (McCabe et al., 2023). The lack of control in teh temporal discretisation (dt), leads to higher numerical errors within the the PRE estimates. fig. 7, we visualise the evaluation of finite difference in 2D+time as 3D con- volution. The finite difference stencil i.e. the convolutional kernel has unit discretisation of da,, dy and dt associated with the problem and is applied output from the neural-PDE spanning the domain æ, y, where r E X), E [0, Y. t E (]. Though the dis- cretisation deployed within CP-PRE stems from the neural PDE solver, consistently delivers guaranteed coverage regardless of resolution. Even with coarser discretisation, CP-PRE remains valuable as it allows for: statistical identi- fication of poorer fit regions (marginal formulation), high- lights physically inconsistent predictions (joint formulation) and enables relative assessment of physical inconsistencies across predictions. While residuals may be inflated with coarser discretisation, the corresponding bounds reflect this inflation, preserving the relative information about physical inconsistency across series of predictions.
0.005
D(u)
0.000
0.005
-0.010
00.015
0.0
0.5
1.0
1.5
2.0
X
(a) CP-PRE over coarser discretisation
Fine (dt=0.005)
the signal i.e. the
PRE
0.010
Upper Bound
0.005
D(u)
0.000
-0.005
0.0
0.5
1.0
1.5
2.0
X
(b) CP-PRE over finer discretisation
Coverage
1.0
ideal
Y-ais
Empirical Coverage
Fine
X-éxia
0.6
0.2
'o.6"
0.8
1-a
(c) Guaranteed coverage irrespective of discretisation emror.
Figure 8. CP-PRE provides guaranteed coverage irrespective the discretisation associated with the model outputs., however, the width the obtained coverage bounds indicates discretisation error associated with the gradient estimation. Coverage taken for a 90% coverage.
Figure 7. PRE estimation using the 3D convolutions with finite difference stencils as convolutional kemels being applied over neural-PDE predictions as the signals.
20

---

Calibrated
UQ
E. Initial and Boundary Conditions
As mentioned Section 4.3, the focus of our experiments has been in quantifying the misalignment of the model with the PDE in the domain of the problem. A well-defined PDE is characterised by the PDE the domain, the initial condition across the domain at / O and the boundary conditions, reflecting the physics at the boundary. Within a neural-PDE setting, the initial condition does not need to be enforced or measured for as the neural-PDE is set up as an initial-value problem, taking in the initial state to au- toregressively evolve the later timesteps and hence does not come under the purview of the neural-PDE's outputs. The boundary conditions, whether Dirichlet, Neumann peri- odic, follows residual structure outlined in Equation (2), allowing us to use it as PRE-like nonconformity score for performing conformal prediction. In all the problems we have under consideration, the PDEs are modelled under periodic boundary conditions:
du OX O; XE 8Q
(26)
By deploying the eqn 26 as the PRE across the boundary, we can obtain error bars over the boundary conditions as well. Within fig. 9, we demonstrate the error bars obtained by using the boundary conditions as the PRE nonconformity scores for the Navier-Stokes equations.
Margina CP
Jo
Coverogge
Marginal
Lower
Ideal Manginal
Figure 9. Error bars obtained the boundary conditions over the right wall of domain of the Navier-Stokes Equation using Marginal and Joint CP. The empirical coverage obtained using the boundary condition as the PRE nonconformity score is also given.
21

---

Calibrated Physics-Informee UQ
F. Toy Problems: 1D cases
it provides us with mechanism to perform rejection criteria for predictions. Within joint-CP, bounds dictating 1 coverage suggest that approximately, x 100% predictions from the same domain will fall outside the bounds and can be rejected. Further details about the physics, parameterisa- tion of the initial conditions, model and its training can be found in Appendix G.
F.1. Advection Equation
Consider the onc-dimensional advection equation
Ou au
+ 0x
O.
(27)
at
The state variable of interest is bounded within the domain I 20, E [0, 0.5] and moves within the domain constant velocity v. Each trajectory is generated by solving Equation (27) using Crank-Nicolson method (Crank & Nicolson, 1947). Data is sampled using a parameterised initial condition that characterises the amplitude and posi- tion of the Gaussian field. Generated data is used to train a ID FNO that takes in the initial condition and autoregres- sively with a step size of 1, learns to map the next 10 time frames outlined in Equation (28). A reproducible script is attached to the supplementary material.
F.2. Burgers Equation
Consider the ID Burgers' Equation
du du pu Ox2'
Ot + u dz = v
(29)
The state variable of interest u is bounded within the domain & E [0, 2), E 1.25. The field is prescribed by a kinematic viscosity v = 0.002. Data is generated by solving Equation (29) using spectral method (Canuto et al., 2007). Data sampled using a parameterised initial condition is used to train a ID FNO that takes in the initial distribution of the state and learns to autoregressively predict the PDE evolution for the next 30 time frames.
n ff"")
(28)
represents the FNO and nn and u = represents the cur- rent (input) and the predicted future state of the system (output).
3.0
0.10
Prodiction
0.20
0.00
PAE
Marginal Narqinal
0.25
Lower join
Upperjeitt
-0.003
-0.05
5.6
-0.002
2.4
8.00
0/3
0.5
1.5
1.5
2.0
5.2
-0.cs
Prodiction
Ground Truhh
Upper Marginal
Lowear Nargital
Lowner Joim Upper (elrt
06
0.5
1/0
3.1
20
2/2
Predictoom
Ground Thuth
Lower Marginal
Upper Mong nol
Lower guit
Upter joine
Figure 11. Burgers' Equation: (Left) Comparing the neural PDE (FNO) performance with that the physics- based numerical solver at the last time instance. (Middle) Upper and lower bounds for 90% coverage obtained by performing marginal-CP. (Right) Upper and lower bounds for 90% coverage guaranteed by joint-CP over the residual space.
Figure 10. Advection Equation: (Left) Comparing the neural PDE (FNO) performance with that of the physics-based numerical solver at the last time instance. (Middle) Upper and lower bounds for 90% coverage obtained by performing marginal-CP. (Right) Upper and lower bounds for 90% coverage obtained by perform- ing joint-CP. Marginal-CP provides tighter bounds for a prediction as opposed to joint-CP, whereas joint-CP provides method of employing a relative sense of reliability of a prediction within a domain.
Figure 11 illustrates the guaranteed bounds over the residual space of Equation (29) using marginal and joint-CP formu- lations for 90% coverage. Marginal-Cr CP provides cell-wise coverage, yielding tighter bounds for specific subdomains. Joint-CP provides bounds 50 times larger than that of the marginal-CP as it covers the entire prediction domain. De- spite the large bounds, approximately a x 100% predictions fall outside it as given in Figure 3. For details on physics, initial condition parameterisation, model, and training, see Appendix H.
Figure 10 demonstrates the guaranteed bounds obtained over the residual space of Equation (27) with both the marginal and Joint-CP formulations. Being cell-wise, marginal-CP guarantees coverage for each discretised point within the spatio-temporal domain. This allows for tighter bounds and error quantification of interested subdomains within regions but does not provide any guarantee across the entire predic- tion domain. Joint-CP acts across the entire domain and provides a guarantee as to whether a prediction (instead of a single cell) will fall within the domain or not. Larger bounds are observed as they extend over the multivariate nature of the output. Though this comes with bigger bounds,
22

---

Calibrated Physics-Informee UQ
G. 1D Advection Equation
Table 7. Architecture of the 1D FNO deployed for modelling ID
G.I. Physics
Advection Equation
Part
Layer
Outpur Shape
Consider the one-dimensional advection equation, parame-
terised by the initial condition:
Input
(50, 200, 1)
Lifting
Linear
(50, 200, 16)
Fourier
Tourceerld/coo/////////EEE
(50. 1. 16. 200)
Fourier
Tourzlerlds/con//////tt//EEE
(50, I. 16. 200)
Fourier
footrier1aa/con///////////
(50, 1, 16, 200)
du du Ot òr
uD
x [0, t 0.,5.
Fourier Fonzzierld/Coo///////////
(50, 1. 16, 200)
Fourier Poorieerb//cmm11///dd////..
(50, 1, 16, 200)
u(z t 0)- (e x)²
(30)
Projection
Fourier 6
Linear
(50, 1, 200, 256)
(50, 1, 16, 200)
Projection 2
Linear
(50, 1, 200. 1)
Here defines the density of the fluid, the spatial coordi- nate, / the temporal coordinate and v the advection speed. initial condition is parameterised by A and X, representing the amplitude and position of a Gaussian distribution. A no-flux boundary condition bounds the system.
The numerical solution for the above equation is built us- method
ing a finite difference solver with implemented in Python. We construct a dataset by perform- ing& Latin hypercube sampling across parameters A, X. Each parameter is sampled from within the domain given in Table 6 to generate 100 simulation points, each with its own initial condition. Each simulation is run for 50-time iterations with a At = 0.01 across a spatial domain span- ning [0,2], uniformly discretised into 200 spatial units in the x-axis.
Figure 12. Advection Equation: Marginal-CC with a 0.5
G.3. Calibration and Validation
To perform the calibration as outlined in Section 5, model predictions are obtained using initial conditions sampled from the domain given in Table The same bounded do- main for the initial condition parameters is used for cali- bration and validation. 100 initial conditions are sampled and fed to the model to obtain and prediction for both the calibration and the validation.
Table 6. Domain range of initial condition parameters for the 1D) advection equation.
Parameter
Domain
Type
Amplitude (A)
50, 200]
Continuous
Position (X)
[0.5, 1.0]
Continuous
G.2. Model and Training
We one-dimensional FNO model the evolution of the convection -ifffutinnn equation. The FNO learns to per- form the mapping from the initial condition to the next time instance, having a step size of 1. The model autoregres- sively learns the evolution of the field up until the 10th time instance. Each Fourier layer has 8 modes and a width 16. The FNO architecture can be found in Table 7. Considering the field values governing the evolution of the advection equation are relatively small, we avoid normalisations. The model is trained for up to 100 epochs using the Adam opti- miser (Kingma & Ba, 2015) with a step-decaying learning rate. The learning rate is initially set to 0.005 and scheduled to decrease by half after every 100 epochs. The model was trained using an L-Plloss (Gopakumar et al., 2024b).
Figure 13. Advection Equation: joint-CP with a = 0.5
23

---

Calibrated Physics-Informed UQ
H. 1D Burgers Equation
of 32. The FNO architecture be found in Table 9. We employ linear range normalisation scheme, placing the field values between -1 and 1. Each model is trained for to 500 epochs using the Adam optimiser (Kingma & Ba, 2015) with step-decaying learning rate. The learning rate is initially set to 0.005 and scheduled to decrease by half after every 100 epochs. The model was trained using an LP-loss (Gopakumar et al., 2024b).
H.1. Physics
Consider the onc-dimensional Burgers' equation:
0u du
8au + dx dr
= 1/
Øt
(x, t = 0) = sin(arT) + cos(-Brx)) + cosh(ytT))
Table 9. Architecture of the ID FNO deployed for modelling ID Burgers' equation
(31)
Part
Layer
Ouipui Shape
where u defines the field variable, y the kinematic viscos- ity, æ the spatial coordinate, the temporal coordinates. 3 and y are variables that parameterise the initial condi- tion of the PDE setup. The system is bounded periodically within the mentioned domain.
Input
(50, 1. 200, 1)
Lifting
Linear
(50, 1, 200, 32)
Fourier
(50, 32, 200)
Fourier 2
Touriee22i/Cam1oo/Att//EEE
(50, 32, 200)
Fourier
Tourcier2t/Cunv/t//aat/// (50, 32, 200)
Fourier Tourreer2d/Conn/l///at///
(50, 1, 32, 200)
The solution for the Burgers' equation is obtained by de- ploying a spectral solver (Canuto et al., 2007). The dataset is built by performing Latin hypercube scan across the defined domain for the parameters a; 3, y, sampled for each simulation. We generate 1000 simulation points, each one with its initial condition and use for training.
Fourier Fourler2a/CCav/l//ttte
(50, 32, 200)
Fourier woorier225/Can..d//ddd//
(50, 1, 32, 200)
Projection Linear
(50, 200, 256)
Projection Linear
(50, 200 I)
H.3. Calibration and Validation
To perform the calibration as outlined in Section 5, model predictions are obtained using initial conditions sampled from the domain given in Table The same bounded do- main for the initial condition parameters is used for calibra- tion and validation. 1000 initial conditions are sampled and fed to the model perform the calibration and 100 samples are gathered for performing the validation.
The physics of the equation, given by the various coefficients is held constant across the dataset generation throughout as given in Equation (31). Each data point, in each simulation is generated with different initial condition as described above. The parameters of the initial conditions are sampled from within the domain as given Table 8. Each simulation is run for 500-time iterations with At 0.0025 across spatial domain spanning [0, 2], uniformly discretised into 1000 spatial units in the and y axes. The temporal domain is subsampled to factor in every 1h time instance, while the spatial domain is downsampled to every 5th instance.
Table 8. Domain range of initial condition parameters for the 1D Burgers' equation.
Parameter
Domain
Type
[3.3]
Continuous
3
3,
Continuous
7
3. 3]
Continuous
H.2. Model and Training
We train ID FNO to map the spatio-temporal evolution of the field variables. We deploy an auto-regressive structure that performs time rollouts allowing us to map the initial distribution recursively up until the thth time instance with a step size of 1. Each Fourier layer has 8 modes and width
24

---

Calibrated Physics-Informet UQ
Marginal CP (75% coverage)
=
t= 10
0.015
Lower Upper
0.0100
Lower Uppar
0.0075
0.005
0.0050
0.000
0.0025
0.0000
-0.010
0.00
0.5
0.50
0.75
1.00
1.25
1.50 1.75
2.00
0.00
0.25
0.50
0,75
1.00
1:25
1.50
1.75 2.00
t=
t=
Residual
Residual
Lower Upper
Upper
0.000
-0.001
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1:75
2.00
0.00
0.25
0.50
0.75
1.00
1:25
1.50
1.75
2.00
Figure 14, Burgers Equation: Marginal-CP with a 0.75
Joint (75% coverage)
t 10
0.3
Lower
0.2
mmm MA M
Lower Upper
Mpm Mn rur a Manla
Upper
0.1
0.0
0.0
-0.1
-0.2
-0.2-
-0.3
0.00
0.28
0.50
0.75
1.00
1.25
1.50
1.75
2.00
a.00
0.25
0.50
0.75
1.00
3.25
3.50
1.75
2.00
20
Residiall Upper
Resictual
MM MMMMMMM.. MM Woahh
Anw And Mw
And
0.100
A. ha
Upper
0.05
My
/ p youp
-0.05
-0.05
0.10
-0.15
0.25 1.25
1.50
1.75
0.50
X
Figure 15. Burgers Equation: joint-CP with a 0.75
25

---

Calibrated Physics-Infone Informeed UQ
I. Utilising CP-PRE as measure of model quality
Performance
While evaluating the performance of neural-PDE, it is im- portant to their fit not just to the data but to the underlying physics. CP-PRE will provide guaranteed coverage irregard- less of the quality of the model. It will have considerably wider error bounds when the neural-PDE (whether PINN or Neural Operator) fails comply with the physics. How- ever, believe that this is advantage of our method. In CP-PRE formulation, the bounds are estimated across the PDE residual, where the ground truth for well-fit solution should always be near zero. If we get wide error bars further away from the 0 for potentially high coverage estimates, it a strong indication that statistically the model violates the physics of interest.
0.8
0.6-
Ground Tht
0.2-
0.0
0.0
0.5
1.0
1.5
2.0
(a) Marginal advection perfor-
mance
PRE-CP over Good Model
PRE Lower Bound
Upper Bound
Consider the example with the Advection equation. We have two models, a well-fit (good model) and poorly fit one (bad model). As shown in fig. 16, though we obtain guaranteed coverage in the case of both the bad and good models, the width of the error bars indicates the quality of the model. Taken for 90 % coverage, the width of the coverage bounds obtained over the bad model is substantially larger than that obtained by the good model.
0.01
D(w)
0.00
-0.02
0.0
0.5
1.0
1.5
2.0
X
(b) Good marginal advection
There still could be a concern as to what width can be con- sidered to be within good range within the residual space. This could be estimated by running the PRE convolution operator(s) across single numerical simulation of the inter- ested physics, thereby estimating the impact of the operator in estimating the residual. The PRE over the simulation data will allow us judge what ranges for the coverage width differentiate between a "bad" and "good" model.
PRE-CP over Bad Model
0.075
PRE
Lower Bound
0.050
Upper Bound
0.025-
D(o)
0,000
0.025-
-0.550-
00075-
0.0
0.5
1.0
1.5
2.0
(c) Bad marginal advection
Coverage
1.0
Empirical Coverage
0.8-
0.6
0.4-
02
0.4
0.6
0a
1-a
(d) Coverage analysis
Figure 16. CP-PRE provides guaranteed coverage irrespective of the model performance, however, the width the obtained cov- erage bounds indicates the accuracy of the model in obcying the underlying physics. Coverage taken for a = 0.1 ~ 90% coverage.
26

---

Calibrated Physics-nnformed UQ
J. 2D Wave Equation
J.2. Model and Training
J.1. Physics
We deploy an auto-regressive FNO that performs time roll- outs allowing
map the initial distribution recursively
Consider the two-dimensional wave equation:
up until the 20th time instance with step size of 1. Each Fourier layer has 16 modes and width of 32. The FNO architecture can be found in Table 10. We employ a linear range normalisation scheme, placing the field values be- tween and Each model is trained for up to 500 epochs (32) using the Adam optimiser (Kingma & Ba, 2015) with step-decaying learning rate. The learning rate initially set Ōu( U. 0) 0, u( 2, E 02, 1 (0t500.00 and scheduled to decrease by half after every 100 epochs. The model was trained using an LP-loss (Gopaku- mar et al., 2024b). The performance of the trained model can be visualised in Figure/ 17.
Ou Pu 0u = =0, 0t2
= 0r
+
0y²
2.U 2,VE(-1,1], t [0,1],
(12,, V/, 0) E A (a- x))-(v-))())
0L
(33)
where u defines the field variable, the wave velocity, æ and y the spatial coordinates, the temporal coordinates. A, and Y variables that parameterise the initial condition
Table 11, Architecture the 2D FNO deployed for modelling the 2D wave equation
of the PDE setup. There exists additional constraint to the PDE setup that initialises the velocity of the wave to O) The system is bounded periodically within the mentioned domain.
Part
Layer
Output Shape
Input
(50, 64, 64, I)
Lifting
Linear
(50, 1, 64, 64, 32)
Fourier
Rourier2d//Cnvzdd/Adee
(50, 32, 64, 64)
The solution for the wave equation obtained by deploy- ing spectral solver that uses leapfrog method for time discretisation and Chebyshev spectral method tensor product grid for spatial discretisation (Gopakumar et al., 2023a). The dataset is built by performing Latin hyper- cube scan across the defined domain for the parameters A, X, Y, which accounts for the amplitude and the location of the 2D Gaussian. We generate 1000 simulation points, each one with its initial condition and use it for training.
Fourier Fourier Fourtee22c/Conv22//dd///
(50, 32, 64, 64)
(50, 32, 64, 64)
Fourier Fouriar726/Com/2///dde
(50, 32, 64, 64)
Fourier Paurieer2o//oov22///d/////
(50, 32, 64,
Fourier 6
(50, 32, 64, 64)
Projection
Linear
(50, 64, 64, 256)
Projection 2
Linear
(50, 1. 64,64 1)
b-10
1-21
The physics of the equation, given by the various coeffi- cients, is held constant across the dataset generation through- out, as given in Equation (32). Each data point, in each simulation generated with different initial condition as described above. The parameters of the initial condi- tions are sampled from within the domain as given in Ta- ble 10. Each simulation is run for 150-time iterations with At 0.00667 across spatial domain spanning 1|² uniformly discretised into 64 spatial units in the x and y axes. The temporal domain is subsampled to factor in every 5th time instance only.
84
Figure 17. Wave Equation: Temporal evolution of field associated with the wave equation modelled using the numerical spectral solver (top of the figure) and that of the FNO (bottom of the figure). The spatial domain is given in Cartesian geometry.
Table 10. Domain range of initial condition parameters for the 2D equation.
J.3. Calibration and Validation
To perform the calibration as outlined in Section 5, model predictions are obtained using initial conditions sampled from the domain given in Table 10. The same bounded domain for the initial condition parameters is used for cali- bration and validation. 1000 initial conditions are sampled and fed to the model to perform the calibration and 100 samples are gathered for performing the validation.
Parameter
Domain
Type
Amplitude (A)
[10, 10)
Continuous
X Position (X)
(0.1, 0.5]
Continuous
Y Position (X)
(0.1, 0.5]
Continuous
27

---

Calibrated Informed UQ
Uval
x10-1
PRE: D(uval)
X10-3
3.0
2.5
2.5
2.0
1.5
2.0
1.0
1.5
0.5
0.0
11.0
-0.5
0.5
1.0
0.0
X
X
(a) Ground Truth
(b) PRE over Ground Truth
U preed
x10-1
PRE: D(Upred)
4.0
x10-3
1.0
3.5
0.5
3.0
0.0
2.5
-0.55
2.0
-1.0
1..5
X
X
(c) Prediction from neural PDE
(d) PRE over Prediction
Figure 18. Analysing the PRE the ground truth and the prediction. Though the neural PDE solver is capable of learning seemingly indistinguishable emulation of the physics while exploring the PRE over each tells different story. As opposed to the smooth Laplacian of the PRE over the ground truth, PRE over the prediction indicates noisy solution, potentially arising due to the stochasticity of the optimisation process.
28

---

Calibrated Physics-nfformed UQ
K. 2D Navier-Stokes Equations
1. Each Fourier layer has modes and width of 16. The FNO architecture be found in Table 13. We employ linear range normalisation scheme, placing the field val- ues between -1 and 1. Each model is trained for up to 500 epochs using the Adam optimiser (Kingma & Ba, 2015) with a step-decaying learning rate. The learning rate ini- tially set to 0.005 and scheduled to decrease by half after every 100 epochs. The model was trained using an LP-loss (Gopakumar et al., 2024b). The performance of the trained model can be visualised in Figure 19.
K.1. Physics
Consider the two-dimensional Navier-Stoker equations:
V V = 0,
o
+ v Vー = ル VP,
with initial conditions:
Table 13, Architecture of the 2D FNO deployed for modelling 2D Navier-Stokes equations
u(a, y, 1 0) - sin(2 u(z, y, 1 0) 3in(473x) x E [-11, 1J,
y 1-1,1, 1},
(34)
Part
Layer
Outpur Shape
(35)
Input
(50, 100, 100,
Lifting
Linear
(50. 1, 100. 100 16)
where u defines the x-component of velocity, defines the y-component of velocity. The Navier-stokes equations solve the flow of an incompressible fluid with a kinematic viscos-
Fourier 1
Tburzlee28/Ccan/2///////
(50, 1, 16, 100, 100)
Fourier 2
Bairie22f/Can/2///d///ll
(50, 16, 100. 100)
Fourier
Eoariiez26/Ccanr2///e///uuuu
(50,
Fourier 4
Fouriee228/conn/2///ee////
(50. 16, 100, 100)
ity v. The system is bounded with periodic boundary condi- tions within the domain. The dataset is built by performing Latin hypercube scan across the defined domain for the pa- rameters a, 3, which parameterises the initial velocity fields for each simulation. We generate 500 simulation points, each one with its initial condition and use it for training. The solver is built using a spectral method outlined in Philip Mocz's code.
Fourier 5
Sburier22d/Conn////dddtt
(50, 1, 16, 100, 100)
Fourier 6
Shurier2dd/Conv2///eet////
(50, 16, 100, 100)
Pmnjecton|
Iihear
(50, 100, 100, 256)
Projectio Linear
(50, 100, 100
K.3. Calibration and Validation
To perform the calibration as outlined in Section 5, model predictions obtained using initial conditions sampled from the domain given in Table 12. The same bounded domain for the initial condition parameters is used for cali- bration and validation. 1000 initial conditions are sampled and fed to the model to perform the calibration and 100 samples are gathered for performing the validation.
Each data point, as each simulation is generated with a different initial condition as described above. The param- eters of the initial conditions are sampled from within the domain as given in Table 12. Each simulation is up until wallclock time reaches 0.3 At = 0.001. The spatial domain is uniformly discretised into 400 spatial units in the and axes. The temporal domain is subsampled to factor in every 1h time instance, and the spatial domain downsampled factor every 4th time instance leading to 100 100 grid for the neural PDE.
Table 12. Domain range of initial condition parameters for the 2D Navier-Stokes equations
Parameter
Domain
Type
Velocity x-axis (Ma)
[0.5, 1.0]
Continuous
Velocity y-axis (vo)
[0.5, 1.0]
Continuous
K.2. Model and Training
We train 2D multivariable FNO to map the -patio-temporal evolution of the field variables (Gopakumar et al., 2024b). We deploy an auto-regressive structure that performs time rollouts allowing us to map the initial distribution recur- sively up until the 20th time instance with a step size of
29

---

Calibrated UQ
t=10
t=21
0.6
0.2
Solution
0.0
-0.2
0.0
-0.4
-0.2
9.50
9.75
10.00
10.25
10.50
0.6
0.4
0.0
0.
-0.2
0.05
0.0
02
-0.4
-0.2
9.50
9.75
10.00
10.25
10.50
Spatio-temporal evolution of the horizontal component of velocity (u)
t=10
t=21
0.6
0.2
0.4
Solution
0.0
0.0
-0.2
-0.5
0.0
-0.4
-02
9.50
9.75
10.00
10.25 10.50
0.50
0.6
0.4
PNO
0.0
0.0
0.00
-.2
-0.5
0.25
-0.4
-0.50
-02
9.50
9.75
10.00
10.25
10.50
(b) Sptio-temporal evolution of the vertical component of velocity (v)
t=10
t=21
0.6
0.50
0.2
0.25
Solution
02
0.0
0.0
0.0
-02
-0.25
-0.4
-04
9.50
9.75
10.00
10.25 10.50
0.50
0.5
0.25
0.4
0.2
PNO
0.0
0.00
-0.2
-0.5
-0.25
-02
-0.4
-0.50
0044
9.50
9.75
10.00 10.25
10.50
(c) Spatio-temporal evolution of the pressure field (P)
Figure 19. Navier-Stokes Equations: Temporal evolution of velocity and pressure modelled using the numerical spectral solver (top of the figure) and that of the FNO (bottom of the figure)
30

---

Calibrated Informeed UQ
Marginal CP (+q)
PRE: D Dom (U, v, P)
x10
x1o-
8
Joint CP (+q x mod)
N
3.0
5
2.5
2.0
0
-4
1.5
1.0
-
0.5
X
X
PRE:
Momentumn Equation
(b) Marginal Bounds
(c) Joint Bounds
Figure 20. Navier-Stoken CP using the Momentum Equation (13) as the PRE for neural PDE surrogate model trained to model fluid dynamics. Figure 20(a) depicts PRE, Figure 20(b) depicts upper error bar, marginal for each cell, while Figure 20(c) indicates the upper error bar obtained across the entire prediction space. Both estimated for 90% coverage.
PRE: D Doont (U, v)
Marginal CP (+q)
Joint CP (+q x mod)
1.0
1.2
x102
1.25
1.0
0.5
1.00
0.0
0.6
0.75
-055
0.4
0.50
1..
0.2
0.25
X
(a) PRE of the Continuity Equation over the (b) Upper error bar indicating 90% coverage (c) Upper error bar indicating 90% coverage FNO prediction
with marginal-CP
with joint-CP
Figure 21. Navier-Stokes: CP using the Continuity Equation (12) the PRE for neural PDE surrogate model trained to model fluid dynamics. Figure 20(a) depicts PRE, Figure 20(b) depicts upper error bar, marginal for each cell, while Figure 20(c) indicates the upper error bar obtained across the entire prediction space. Both estimated for 90% coverage.
31

---

Calibrated
UQ
L. 2D Magmetolydrodynamies
Consider the Ideal MHD equations 2D:
sin(2atY)
(36)
V sin(2bt X),
(37)
P 4Cx
(38)
where the density (p), velocity field (v [u, e|) and the pressure of plasma is modelled under a magnetic field
(B [z. Byl) across spatio-temporal domain 2, y [0. LE [0, 5]. 140 is taken to be the magnetic perme- ability of free space. The system is bounded with periodic boundary conditions within the domain. The dataset is built by performing Latin hypercube scan across the defined domain for the parameters a, b, c, which parameterises the initial velocity fields for each simulation. We generate 500 simulation points, each one with its initial condition and use it for training. The solver is built using finite volume method outlined in Philip Mocz's code.
Each data point, as in each simulation is generated with a different initial condition as described above. The param- eters of the initial conditions are sampled from within the domain given in Table 12. Each simulation is run up until wallclock time reaches 0.5 with varying temporal discreti- sation. The spatial domain is uniformly discretised into 128 spatial units in the x and y axes. The temporal domain is downsampled to factor in every 25th time instance.
0p + 4. (0 = 0,
A9
B x (v B) VP.
10
140
d
P
dt
p7
Table 14, Domain range of initial condition parameters for the 2D MHD equations
Ø
= A x (V B),
0t
Parameter
Domain
Type
B B= 0,
Velocity x-axis (a)
[0.5, 1.0
Continuous
Velocity y-axis (6)
[0.5, 1.0
Continuous
Pressure (c)
[0.5, 1.0
Continuous
L.1. Model and Training
We train 2D multi-variable FNO map the spatio- temporal evolution of the 6 field variables collectively, We deploy an auto-regressive structure that performs time roll- outs allowing us to map the initial distribution recursively up until the 20th time instance with a step size of 1. Each Fourier layer has modes and width 16. The FNO architecture can be found in Table 15. We employ linear range normalisation scheme, placing the field values be- tween -1 and Each model is trained for up to 500 epochs using the Adam optimiser (Kingma & Ba, 2015) with step-decaying learning rate. The learning rate is initially set to 0.005 and scheduled to decrease by half after every 100 epochs. The model was trained using an LP-loss (Gopaku- mar et al., 2024b). The performance of the trained model
with initial conditions:
32

---

Calibrated Informed UQ
be visualised in Figures 22 and 23.
Table 15, Architecture of the 2D FNO deployed for modelling 2D
MHDD equations
Part
Layer
Output Shape
Inpot
(50, 128, 1)
Lifing
Linear
(50, 128, 128 16)
Fourier
TourLer2d//oov22///4/// GELU
(50, 1. 16, 128, 128)
Fourier Pourier2dd/Cav2//A44///E//
(50, 16, 128, 128)
62
Fourier Vourier2d/Conv2///dd GEUU
(50, 1, 128, 128)
Fourier Vourier2dd/Co//////6///
(50,
128)
Fourier Vourier2dd/Con/2///kd//
(50, 1. 16, 128, 128)
-05
Fourier
Pourier2dd/Cow/////dd////
(50, 1, 16, 128, 128)
Projection
Linear
(50, 1, 128, 128, 256)
Projection
Linear
(50, 128
-05
L.2. Calibration and Validation
To perform the calibration as outlined in Section 5, model predictions are obtained using initial conditions sampled from the domain given in Table 12. The same bounded domain for the initial condition parameters is used for cali- bration and validation. 100 initial conditions are sampled and fed to the model to perform the calibration and 100 samples are gathered for validation.
(a) Spatio-temporal evolution of density (p)
-05
M0000
02
05
2000
(b) nentioltemporal evolution of the horizontal component of veloc- ity (u)
52
-64
68
.2
23
-04
-06
-03
Spatio-temporal evolution of the vertical component of velocity
Figure 22. MHD Equations: Temporal evolution of velocity and
pressure modelled using the numerical solver (top of the figure) and that of the FNO (bottom of the figure). (Continued on next page)
33

---

Calibrated Inforreed UQ
676
820
626
04
821
IICE
(a) Spatio-temporal evolution of the pressure field (P)
64-
30
675 1000 WUS
30
1000
(b) tentora evolution of the horizontal component of the magnetic field (B=)
806
8.0
10
+0.00
-41
-22
10.00
62
902
AL
-23
8.00
30
I0CE
(c) Spatio-temporal evolution of the vertical component of the mag- netic field (By)
Figure 23. MHD Equations: Temporal evolution of velocity and pressure modelled using the numerical solver (top of the figure) and that of the FNO (bottom of the figure). (Continued from previous page)
34

---

Calibrated Physics-Informit UQ
Marginal
joint
Marginal
0.04
PA Lower
Lower Upper
Upper
6.62
0..5
PRS: Doa.sei1...
0.00
0.00
0.00
-0.02
-b2
-0.0
Induction Equation
Energy Equation
Figure 24. MHD: Slice plots along the x-axis (sliced at = 0.5m) indicating the marginal and joint coverage (90%) obtained over the neural PDE modelling MHD equations using the induction equation Equation (17) (on the left) and energy equation Equation (16) (on the right). Marginal coverage, evaluated cell-wise, generates tight bounds the PRE, whereas joint coverage spanning across the spatio-temporal domain introduces wider bounds.
PRE: DinductionV,, B)
Marginal CP (+q) x10-1
Joint CP (+ x mod)
x10-2 4
1.25
6
2
1.00
4
0
0.75
2
0.50
2
0.25
x
X
(a) PRE of the Induction Equation Equa-(b) Upper error bar indicating 90% cover-(c) Upper error bar indicating 90% coverage tion (17) over the FNO prediction
age with marginal-CP
with joint-CP
Figure 25. MHD:: CP using the Induction Equation (17) as the PRE for neural PDE surrogate model solving the Ideal MHD equations. The time instance the prediction shown.
PRE: D epn V.P, B) x
Marginal CP (+q) x10
Joint CP (+ q mod)
1.0
2.0
0.8
5
1.5
0.6
1.0
00
0.4
-0.5
0.2
-5
X
X
X
(a) PRE of the Energy Equation Equa-(b) Upper crror bar indicating 90% coverage (c) Upper error bar indicating 90% coverage with joint-CP
tion (16) over the FNO prediction
with marginal-CP
Figure 26. MHD: CP using the Energy Equation (16) as the PRE for neural PDE surrogate model solving the Ideal MHD equations. The last time instance of the prediction is shown. Error bars obtained using joint CP are an order of magnitude higher than that obtained by marginal CP as it is measured across the entire spatio-temporal domain rather than for each cell.
35

---

Calibrated Inforreed UQ
Marginal CP (+q)
Joint CP (+q x mod)
PRE: Dcant. (p. V)
x10
100=l
-1.5
1.5
2.0
1.0
1.0
1.5
0.5
0.0
0.5
1.0
-0.5
0.5
x
X
(a) PRE of the Continuity Equation Equa-(b) (b) Upper error bar indicating 90% cover-(c) Upper enror bar indicating 90% cover- tion (14) over the FNO prediction
age with marginal-CP
age with joint-CP
Figure 27. MHD: CP using the Continuity Equation (14) as the PRE for neural PDE surrogate model solving the Ideal MHD equations.
Marginal CP (+q)
PRE: D Dausse
Joint CP (+q x mod) x
6
0
1.0
2
00.5
PRE of the Divergence Equation Equa-(b) Upper bar indicating 90% cover-(c) Upper bar indicating 90% age with joint-CP
tion (18) over the FNO prediction
age with marginal-CP
Figure 28. MHD: CP using the Gauss's law for magnetism Equation (18) the PRE neural PDE surrogate model solving the Ideal MHD equations.
36

---

Calibrated Physics-Informun UQ
M. Plasma Modelling within a Tokamak
M.1. Physics and Problem Setting
Ideal Residual
We evaluate our uncertainty quantification framework on simplified Magneto-hydrne hydrodynamies (MHD) model in toroidal geometry, leveraging the dataset and neural archi- tecture from (Gopakumar al., 2023b). The dataset is generated using the JOREK code (Hoelzl et al., 2021) with
0.8
Empirical Coverage
physics model similar to the Reduced-MHD model but with electrostatic and isothermal constraints. In this setup, only the density p. electric potential +, and temperature T fields are evolved.
0.6
0.4
The physical system models the dynamics of toroidally axisymmetric density blob initialized on top of low back- ground density. Without plasma current maintain confine- ment, the pressure gradient in the momentum equation acts as a buoyancy force, causing radial blob motion. This sim- plified scenario serves as proxy for studying the evolution of plasma filaments and LocelLoo Modes (EI.Ms) in tokamak devices.
0.2
0.2
0.4
0.6
0.8
1-alpha
Figure 29. Empirical coverage obtained by performing CP-PRE over JOREK-PNO
M.2. Dataset Generation
The dataset consists of 120 simulations (100 training, testing) generated by solving the reduced MHD equations using JOREK with periodic boundary conditions. The ini- tial conditions vary in the blob's position, width, and am- plitude. Each simulation is performed within the toroidal domain with 1000 timesteps, downsampled to 100 slices and 200x200 bi-cubic finite-element spatial grid, downsampled 100
predictions orders magnitude faster than the numerical solver. For a complete description of the MHD system and additional experimental details, we refer readers to (Gopaku- et al., 2023b).
M.4. Calibration and Validation
For this system, we compute physics residual errors (PRE) for the temperature equation (equation from (Gopakumar et al., 2024b)), it comprises all the variables modelled by the neural PDE. Empirical coverage obtained by performing CP-PRE over the JOREK FNO is shown in Figure 29.
M.3. Model Architecture and Training
We employ the Fourier Neural Operator (FNO) architec- ture from (Gopakumar et al., 2023b) with the following specifications:
Input: 20 time instances of field values 100x100 grid
Autoregressive prediction up to 70 timesteps
Output: Next time instances
Fourier layers with width 32 and 16 modes
Physics Normalisation followed by standard linear range normalization to [-1,1]
The training was conducted on single A GPU with the Adam optimizer having an initial learning rate of 0.001 with halving every 100 epochs for 500 epochs total using a relative LP loss function. This model achieves normalised
MSE of 4e on each field variable while providing
37

---

Calibrated Physics- Informedd UQ
N. Magnetic Equilibrium in a Tokamak
P1
P6
10
N.1. Physies and Problem Setting
tokamak uses magnetic fields to confine and shape a plasma for nuclear fusion. While the main toroidal field (TF) running the long way around the torus provides the primary confinement, the poloidal field (PF) coils running in loops around the cross-section are crucial for plasma stabil- ity and performance. These coils serve multiple functions: they induce the plasma current that generates an additional poloidal field component, shape the plasma cross-section, and provide vertical stability control. Without these care- fully controlled poloidal fields, the plasma would be un- stable and quickly lose confinement. The structure of the tokamak with emphasis on the magnetic fields and the coils are given in Figure 31.
+S1
P9
6
-22
-P3
0
The shape the plasma cross-section significantly impacts performance, with modern tokamaks typically using a "D- shaped" plasma with strong elongation and triangularity. This shape, controlled actively by varying currents differ- ent poloidal field coils, is superior simple circular cross- section as it allows for higher plasma current at given magnetic field strength, improving confinement. The D- shape with triangularity also provides better stability against plasma instabilities, particularly edge localized modes, and enables better access to high-confinement operation. Main- taining this shape requires sophisticated feedback control systems since the plasma shape naturally unstable and needs continuous adjustment, especially for maintaining the separatrix the magnetic surface that defines the plasma boundary and creates the X-point for the divertor. A sample magnetic equilibrium across the poloidal cross-section is showcased Figure 30, with the equilibrium shown in the contour plots and the poloidal field coils indicated with the grey block. The structure of the tokamak and the separatrix are indicated in black and red respectively.
-10
2
R [m]
Figure 30. Sample Equilibrium plot showcasing the magnetic equi- librium as contour plots observed given poloidal field (PF) configuration (PF coil locations are indicated in blue blocks), The poloidal cross-section of the tokamak shown here, with the structural boundary in black and the closed flux surfaces in pink (Meyer Team, 2024). Our problem looks at mapping the equilibrium across the spatial domain for given PF coil location.
Considering the impact of the PF coils on the equilibrium and, hence on the plasma efficiency, its placement within the structure is an open design problem for future tokamaks. Traditionally, the design space is explored by obtaining the equilibrium associated with certain PF coil configuration, which involves solving the Grad-Shafranov equation using numerical solvers, rendering them computationally expen- sive. As alternative, we construct surrogate model to explore the design space significantly faster.
while being conditioned on the locations of the PF coils in the latent space. The performance of the model can be found Figure 33, The model was trained on 512 simula- tions of the Grad- Grad-Shafranov simulations obtained using the FreeGSNKE code (Amorisco et al., 2024). The coil currents are kept constant we are looking to identify the ideal coil configuration for a steady-state plasma.
N.2. Surrogate Model
N.3. Calibration and Validation
A conditional auto-encoder, given in Figure 32, is trained as a surrogate model that can approximate the magnetic equilibrium for a given configuration of PF coils. The auto-encoder takes in as input the spatial domain the tokamak the input and outputs the magnetic equilibrium
CP-PRE is performed by evaluating the Grad-Shafranov equation over the solution space of the surrogate model, with inputs generated by sampling within the bounds of the PF coil locations. By using the GS equation, we can identify the predictions from the surrogate model that generate untenable
38

---

Calibrated Physics-Informec UQ
Inner poloidal field coils (primary transformer circuit)
Poloidal magnetic field
Outer poloidal coils
(for plasma positioning and shaping)
Toroidal field coils
Resulting belical magnetic field
Toroidal magnetic field
Plasma electric current (secondary transformer circuit)
Figure 31. A simple schematic diagram of a generic tokamak with all the main magnetic components and fields shown (Li et al., 2014).
The poloidal field coil magnets (grey) are that which this work aims to optimise.
Encoder
Decoder
Ground Truth
Prediction
5
100
10
15
15
R,Z Mesh
Equilbirum
-2
-4
25
25
-6-
-30
30
PF Coil Locations
8-
2
Figure 32. Conditional auto-encoder developed as a surrogate model for mapping the poloidal field coil locations to the responding magnetic equilibrium under constant coil currents.
Figure 33, Comparing ground truth with prediction the surrogate model.
equilibria. This allows to explore the design space quickly while adding trustworthiness to your surrogate model. Both marginal and coverage obtained using CP-PRE with GS equations over the surrogate model is indicated in Figure 34. Marginal-CP shows smooth coverage as it represents the coverage averaged across each cell.
Coverage
Ideal
Marginal
Empirical Coverage
0.8-
Joint
0.4-
o.o
1-d
Figure 34. Marginal and joint empirical coverage obtained by per- forming CP-PRE over the Grad- Shafranov equation
39


Reviews:
[{"review_id": "NDlm5332Bk", "reviewer": "Reviewer_MQ8g", "summary": "Papers consider uncertainty quantification of PDEs. They claim that by utilising a physics based approach they can quantify and calibrate the model’s inconsistencies with the PDE.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# The big picture:\n\nWe appreciate this opportunity to clarify our motivation. Our work stems from the need to make neural PDE solvers more practical for scientific modelling. Numerical PDE solvers have been essential to scientific modelling since the 1950s, enabling cost-effective simulation of complex scenarios. However, they still impose significant computational burdens in complex settings.\n\nNeural PDEs offer a promising middle ground—networks that learn physics approximations could quickly explore design spaces before running expensive numerical simulations. For this workflow to be practical (as shown in Figure 1), we must quantify neural PDE model performance with reliable error bounds.\n\nWhile various Bayesian methods can be applied to neural PDE solvers, they fail to guarantee coverage in high-dimensional settings. Extensions of inductive conformal prediction to spatio-temporal domains require substantial calibration data for error bounds.\n\nOur approach addresses these limitations by using PDE residuals as nonconformity scores, providing guaranteed bounds across conservative PDE variables without requiring data or sampling, functioning post-hoc without modifications. Table 2 in Appendix C qualitatively compares our method's advantages against Bayesian UQ and standard CP methods.\n\n***\n***\n# Utility of this method and the meaning of these figures:\n\nOur approach begins with predictions of spatio-temporal physical field variables from a neural PDE solver. We then transform these predictions to their conservative form using the differential operators that characterise the PDE. This transformation captures how well the solution satisfies conservation equations, identifying regions of physical inconsistency. The Physical Residual Error (PRE) estimates reveal regions where the predicted dynamics deviate from the ground truth. By applying Conformal Prediction (CP) to these PRE values, we obtain calibrated bounds across conservative variables that provide meaningful information to domain experts. The figures illustrate how conservative variables are modelled across space and time and the regions where the model struggles to understand and learn the known physics. \n\nAdditionally, our PRE-CP framework serves as an indicator of model quality, as demonstrated in Appendix I. While PRE-CP guarantees coverage regardless of model fit, the width of the error bars effectively measures model quality. Figure 18 compares (a) predictions from poorly-fit (bad) and well-fit (good) models, showing that (b) error bounds for the well-fit model are substantially narrower than (c) those of the poorly-fit model.\n\n***\n***\n\n# Procedure:\n\nWe thank the reviewer for this helpful suggestion. Due to page constraints, we limited our discussion of conformal prediction in the main text while referencing relevant literature for broader context. We have outlined an abridged algorithmic procedure below and have added a detailed one to the appendix. \n\n1.\tNeural PDE Solver Setup: Define PDE simulator, generate data and train neural approximator\n2.\tPhysics Residual Error: Calculate $|D(\\hat{f}(X))|$ using differential operators as convolutions. \n3.\tCalibration:\n\n  a. Sample ${X_1,...,X_n}$, get predictions ${\\hat{f}(X_1),...,\\hat{f}(X_n)}$\n\n  b.\tCompute PRE scores ${|D(\\hat{f}(X_1))|,...,|D(\\hat{f}(X_n))|}$\n\n  c.\tFor confidence $1-\\alpha$, find $\\hat{q}^\\alpha$ as the $\\lceil(n+1)(1-\\alpha)\\rceil/n$ quantile from the PRE scores. \n4.\tApplication:\n\n  a. For new $X_{n+1}$, calculate $|D(\\hat{f}(X_{n+1}))|$\n\n  b.\tValid if within $\\mathbb{C}^\\alpha = [-\\hat{q}^\\alpha, \\hat{q}^\\alpha]$\n\n5.\tInterpretation: Apply cell-wise for local bounds or use supremum for global bounds\n\nWe appreciate your feedback and thank you for making the paper accessible to a broader audience.\n***\n***\n\n# Comparison to other methods:\n\nOur method's coverage performance is documented in Tables 3-5 in Appendix C, where we compare against both Bayesian UQ approaches (deep ensembles, MC dropout, stochastic-weighted averaging, Bayesian neural networks) and data-driven conformal prediction using absolute error residuals (CP-AER).\n\nOur method (CP-PRE) consistently provides guaranteed coverage across all experiments in both in-distribution and out-of-distribution evaluations, with performance comparable to CP-AER. However, it's important to note that CP-AER requires extensive simulation data for calibration, while our method does not. This difference is quantified in the tables, showing that the computational cost of evaluating PRE is only a fraction of the time needed to acquire new simulation data, but has been explicitly mentioned in response to reviewer X8Gm's questions. \n\nFollowing reviewers' suggestions, we have moved Tables 3-5 to the main paper and hope that demonstrates the advantages of our method. \n\nWe apologise for the confusion that this may have caused and welcome any further clarification requests.", "claims_and_evidence": "(see Other Strengths And Weaknesses))", "methods_and_evaluation_criteria": "(see Other Strengths And Weaknesses))", "theoretical_claims": "(see Other Strengths And Weaknesses))", "experimental_designs_or_analyses": "(see Other Strengths And Weaknesses))", "relation_to_broader_scientific_literature": "PDEs are important in physics.", "essential_references_not_discussed": ".", "other_strengths_and_weaknesses": "They claim that by utilising a physics based approach they can quantify and calibrate the model’s inconsistencies with the PDE. I can follow the paper in high level, meaning (based on my understanding) they calculate score function PRE using the stencil approximation and calculate some kind of conformal prediction (CP). But I have to admit that I do not understand the pig picture meaning that why they are doing that or how these CPs can be used in practice?\n\nMaybe the reason is that I am not previously familiar with CP and cannot completely follow the idea based their description. Probable due to size constraints this CP section is quite short and mostly summarises some key definitions and results from previous works, but it does not explain so much that how the quantities are actually calculated or how those are used. They could explain more that how the approach is actually used. \n\nAlso, results shows PRE (which is basically just the derivate of the solution) and derived CPs for different PDE problems. Those are nice looking pictures, but again overall meaning or how those could help in practice is not clear to me. Furthermore, they also do not compare to any of previous methods.", "other_comments_or_suggestions": "Overall approach or \"big picture\" could be explained better. Perhaps showing an example which shows that how CP is calculated and then how it is used to solve practical problem? Some kind of algorithmic description of the overall approach could also be helpful.", "questions_for_authors": "No further questions", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PI0XSg1myk", "reviewer": "Reviewer_Zd7y", "summary": "This paper proposes a model-agnostic, physics-informed conformal prediction network that provides guaranteed uncertainty estimates independent of input data.", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Quantitative evaluation to baselines\n\nThank you for raising this important point. We'd like to clarify that our framework indeed demonstrates superior performance in guaranteed coverage compared to baseline methods. In Appendix C (Tables 3-5), we comprehensively compare our method (CP-PRE) against standard Bayesian approaches (BNN, MC Dropout, Deep Ensembles, SWA-G) and data-driven inductive conformal prediction using absolute error residual (CP-AER).\n\nOur results show that CP-PRE consistently provides guaranteed coverage across all experiments, including both in-distribution and out-of-distribution evaluations, with performance comparable to CP-AER. While CP-AER achieves similar coverage, it requires extensive simulation data for calibration—a significant limitation.\nFollowing Reviewer X8Gm's suggestion, we've updated our evaluation metrics to include data generation time for CP-AER and residual estimation time for CP-PRE in the revised manuscript, providing a more complete comparison of computational requirements and is shown in the rebuttal to X8Gm. \n\nWe appreciate your feedback and have moved Tables 2-5 to the main text, highlighting methods that achieve estimated coverage for improved clarity. Thank you for helping us enhance the paper's readability.\n\n***\n***\n\n# Missing symbols and improving writing\n\nThank you for your careful review of Section 5. We've thoroughly re-examined this section and have ensured all symbols are properly defined. While Reviewer X8Gm found the writing quality satisfactory, we recognise the importance of clarity for all readers. In our revised manuscript, we've added additional clarification where needed and performed a comprehensive review of all notations to ensure consistency and accessibility throughout the paper. We appreciate your feedback, as it has helped us improve the overall quality of our presentation.\n\n\n\n***\n***\n\n# Aleatoric vs Epistemic Uncertainty\n\nThank you for raising this important conceptual question. The uncertainty quantified by PRE-CP aligns with conformal prediction's characterization of predictive uncertainty.\nFrom one perspective, this could be viewed as aleatoric uncertainty since we construct confidence intervals relative to a specific probability distribution (the distribution from which calibration data, i.e. initial conditions are sampled). Alternatively, it could be considered epistemic uncertainty since we model the neural network's error through confidence intervals (typically used for unknown but fixed quantities).\nWhile we believe the latter interpretation is more appropriate, we acknowledge that the traditional aleatoric/epistemic dichotomy may not be directly applicable to our framework. This distinction is most valuable when both uncertainty types coexist and require separate treatment [1]. Although we could elaborate on this in our manuscript, we felt it somewhat peripheral to our paper's focus on practical uncertainty quantification rather than fundamental statistical theory.\n\n\n[1] S. Ferson, and L. Ginzburg, \"Different methods are needed to propagate ignorance and variability\", Reliability Engineering & System Safety, Elsevier, 1996, https://doi.org/10.1016/S0951-8320(96)00071-3\n\n\n***\n***\n\n\n# Validation plots for all experiments\n\nThank you for this question. We've included complete validation plots and quantitative analysis for all experiments, including plasma modeling and magnetic equilibrium, in the appendices due to space constraints in the main text:\n\n\n•\tPlasma modeling coverage plots appear in Figure 29, Appendix M\n\n•\tMagnetic equilibrium coverage plots are shown in Figure 34, Appendix N\n\n\nWe initially focused on presenting coverage guarantees for widely recognized test cases (Wave, Navier-Stokes, and MHD equations) where we do ablation studies in the main text. Based on your feedback, we will improve cross-referencing to these appendices to ensure readers can easily locate the complete analysis for all experiments.", "claims_and_evidence": "1. The proposed approach is model-agnostic and physics-informed. The physics-informed aspect is evidenced in Section 4, but whether the approach is truly model-agnostic is not clearly demonstrated.\n2. The proposed approach guarantees coverage bounds both marginally and jointly. This is evidenced by Theorem A.1.", "methods_and_evaluation_criteria": "The authors utilize physicist residual error as the nonconformity score, enabling data-free prediction. The proposed approach is validated through comparisons of estimated uncertainty with PDE residuals.", "theoretical_claims": "The authors claim that the proposed approach guarantees coverage bounds both marginally and jointly, as formalized in Theorem A.1.", "experimental_designs_or_analyses": "The authors evaluate the proposed approach using:\n(1) Wave equations,\n(2) Navier-Stokes equation,\n(3) Magnetohydrodynamics,\n(4) Plasma modeling within a tokamak,\n(5) Magnetic equilibrium in a tokamak.\n\nThe estimated uncertainty of the model is compared with the PDE residual.", "relation_to_broader_scientific_literature": "The study contributes to uncertainty quantification (UQ) in physics-informed machine learning and conformal prediction. The proposed method aligns with ongoing research in Physics-informed neural networks (PINNs), neural PDE solvers, and uncertainty quantification in computational physics.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Strengths:\n- The paper introduces an interesting framework that extends conformal prediction to physics-informed models.\n\nWeaknesses:\n1. From the quantitative evaluation in the supplementary material, the proposed approach does not achieve the best performance in most scenarios compared to other baselines. This raises concerns about its effectiveness.\n2. The writing of the paper needs improvement, as many symbols in Section 5 are not clearly defined.", "other_comments_or_suggestions": "- Please highlight the best-performing approach in the quantitative tables in the supplementary material.\n- It may be beneficial to move these tables to the main paper, as quantitative results are important for evaluating the proposed approach.", "questions_for_authors": "1. Uncertainty can be broadly categorized into aleatoric and epistemic. Given that the proposed approach estimates uncertainty independent of input, how does it differ from epistemic uncertainty estimation approaches?\n2. Why are validation plots and quantitative analysis provided only for the first three experiments but not for plasma modeling and magnetic equilibrium?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "pIlfd4rSIY", "reviewer": "Reviewer_Qd74", "summary": "This paper presents a method for estimating uncertainties in neural PDE solvers without requiring labeled data. The authors propose PRE-CP, which combines PDE residuals and conformal prediction. By using the PDE’s own equations as the reference, the method calibrates each model’s physical errors directly. They show that PRE-CP works with standard neural PDE setups (such as wave, Navier-Stokes, and magnetohydrodynamics) and also demonstrate applications to fusion research. Their results indicate that PRE-CP can uncover locations in the model’s predictions where the solution fails to follow the underlying physics, offering a statistically valid way to determine when to trust or question the model’s output.\n\n## update after rebuttal\nRebuttal acknowledged.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# README for Code:\n\nThank you for pointing out this error that happened during anonymisation and giving us a chance to rectify it. For the purpose of the review, we are providing an abridged README below. \n\n## Installation\n\n```bash\npip install -r requirements.txt\n```\n\n## Quick Start\n\nRun standalone experiments (no data or pre-trained models needed):\n\n```bash\n# Marginal bounds for 1D advection\npython -m Marginal.Advection_Residuals_CP\n\n# Joint bounds for 1D advection\npython -m Joint.Advection_Residuals_CP\n```\n\n## PRE Estimation Example\n\n```python\nfrom ConvOps_2d import ConvOperator\n\n# Define operators for PDE: ∂u/∂t - α(∂²u/∂x² + ∂²u/∂y²) + βu = 0\nD_t = ConvOperator(domain='t', order=1)  # time-derivative\nD_xx_yy = ConvOperator(domain=('x','y'), order=2)  # Laplacian\nD_identity = ConvOperator()  # Identity Operator\n\n# Combine operators with coefficients\nalpha, beta = 1.0, 0.5\nD = ConvOperator()\nD.kernel = D_t.kernel - alpha * D_xx_yy.kernel - beta * D_identity.kernel\n\n# Estimate PRE from model predictions\nPRE = D(model(X))\n```\n\n## Advanced Experiments\n\nFor Navier-Stokes, MHD, and other experiments:\n1. Generate data: Run scripts in `Neural_PDE/Numerical_Solvers/`\n2. Train models: Use scripts in `Neural_PDE/Expts/`\n3. Run uncertainty estimation: See scripts in `Marginal/` and `Joint/`\n\n## Repository Structure\n\n- `Joint/`, `Marginal/`: Conformal prediction implementations\n- `Neural_PDE/`: Neural PDE solver implementations\n- `Utils/`: Utility functions\n- `Other_UQ/`: Bayesian Deep Learning benchmarks\n\n\n***\n***\n# Discussion on discretisation:\n\nWe've addressed the impact of spatial/temporal resolution on residual estimates in Appendix D.1, but agree this merits further discussion. The discretisation in PRE-CP stems from the neural-PDE solver itself. With pre-trained models, we typically have limited control over this aspect. The convolutional kernels (and corresponding finite difference stencils) adopt the discretisation present in the predicted data. As shown in Figure 10 for the 1D Advection equation, coarser discretisation leads to inflated error bounds compared to finer discretisation implementations. However, PRE-CP consistently delivers guaranteed coverage regardless of resolution.\nEven with coarser discretisation, PRE-CP remains valuable by:\n\n1.\tProviding statistical identification of poorer fit regions (marginal formulation)\n\n2.\tHighlighting physically inconsistent predictions (joint formulation)\n\n3.\tEnabling relative assessment of physical inconsistencies across predictions\n\nWhile residuals may be inflated with coarser discretisation, the corresponding bounds reflect this inflation, preserving the relative information about physical inconsistency across a series of predictions. We're currently leveraging this property to develop an active learning pipeline for neural PDEs, using PRE-CP formulation as an acquisition function.\n\n***\n***\n# Unknown PDE terms:\n\n\nThank you for this insightful question. Our method fundamentally relies on formulating the nonconformity score as an equality constraint through the residual formulation, as demonstrated in the Theorem in Appendix A. This approach allows us to perform conformal prediction without requiring data.\n\nWhen certain PDE terms are unknown, the equality constraint may not be fully satisfied, potentially limiting PRE-CP's direct applicability. However, PRE-CP doesn't necessarily require complete knowledge of the entire PDE family. As shown in our Navier-Stokes and MHD equation examples, we can still derive meaningful bounds using only one conservation law (continuity or momentum) without explicitly incorporating all equations in the family.\n\nIt's worth noting that our work primarily targets neural PDE surrogates for forward modelling where the PDE terms are known. This focus enables us to provide guaranteed uncertainty quantification for these specific applications.\n\n***\n***\n# Computational Overheads:\n\nThe computational overhead for PDE residual evaluation is relatively low due to our use of highly optimised convolutional operations. The cost hierarchy in our workflow is:\n\n1.\tRunning full simulations (e.g., 350 core hours for Tokamak plasma modelling)\n\n2.\tTraining neural PDE models (e.g., 6 hours on a single A100 for FNO)\n\n3.\tModel inference (e.g., 90 seconds on a standard laptop)\n\n4.\tResidual estimation (less than 5 seconds)\n\n\nWe've optimised our ConvOps library to minimise computational costs through:\n\n\n1.\tAdditive kernels for linear PDE components, reducing the number of required convolutions (detailed in Appendix D)\n\n2.\tSupport for both spectral and direct convolutions, providing speedups for high-resolution grids [1]\n\n\nThese optimisations maintain the physics residual equality constraint and prediction exchangeability, preserving our coverage guarantees.\n\n[1] Rippel, O., Snoek, J., & Adams, R. P. (2015). Spectral Representations for Convolutional Neural Networks. arXiv preprint arXiv:1506.03767. Retrieved from https://arxiv.org/abs/1506.03767", "claims_and_evidence": "The key claim is that using a physics-informed residual as a nonconformity score leads to valid coverage guarantees. The paper supports this with theoretical arguments and empirical validation on multiple PDEs.", "methods_and_evaluation_criteria": "The authors employ standard PDE examples (wave, Navier-Stokes, MHD) and then apply the approach to practical fusion applications. The authors use recognized PDE solvers for reference. They also describe how to generate PDE residuals efficiently via finite difference stencils", "theoretical_claims": "The proof sketches in the appendices appear coherent. A deeper check would require verifying each step of their derivation regarding exchangeability assumptions for PDE initial conditions, but no immediate flaws are evident.", "experimental_designs_or_analyses": "The experiments compare coverage results across wave, Navier-Stokes, and magnetohydrodynamics PDEs, and also show real-world scenarios. The analyses use marginal and joint coverage and show how coverage is empirically measured. The experimental protocols appear logically consistent, with error bars and coverage curves displayed clearly.", "relation_to_broader_scientific_literature": "This work expands conformal prediction into PDE-solving contexts by incorporating PDE residuals. It aligns with recent research on physics-informed methods for PDE-based modeling (e.g., PINNs, neural operators). Unlike many UQ approaches requiring labels or Bayesian sampling, this paper uses physics-based residuals, filling a gap in data-free calibration methods. It also relates to prior CP frameworks for high-dimensional data, extending them to spatio-temporal PDE outputs.", "essential_references_not_discussed": "Most relevant PDE-based operator learning work is cited (Fourier Neural Operators, Physics-Informed Neural Networks), and the paper mentions conformal prediction references for spatio-temporal data.", "other_strengths_and_weaknesses": "Strengths:\n+ It provides both marginal and joint coverage formulations, giving users the option to pinpoint local errors or reject predictions at a global scale.\n+ It provides coverage guarantees without labeled data.\n+ The appendices are thorough, showing proofs for the theoretical aspects and comparisons across multiple UQ methods and PDE scenarios.\n\nWeaknesses：\n+ The approach requires knowing the PDE precisely, so it is not suitable if the physical model is partially unknown or has uncertain terms.\n+ The residual-based method can be sensitive to discretization, especially in the temporal domain, potentially inflating error bars in coarser meshes.\n+ Initial and boundary conditions, while addressed, could receive a clearer explanation in the main text for completeness.", "other_comments_or_suggestions": "+ README for Code: Including a short README in the supplementary code would help others replicate the experiments.\n+ Extended Discretization Analysis: More discussion about how changing resolution in space/time affects the residual estimates would be helpful.", "questions_for_authors": "+ If certain PDE terms or coefficients are unknown, could this approach handle partially specified physics by defining an approximate residual operator?\n+ How large is the overhead when evaluating PDE residuals for high-resolution spatio-temporal domains, and are there workarounds to reduce it without losing coverage guarantees?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "uoaUX4XMa7", "reviewer": "Reviewer_X8Gm", "summary": "The paper focuses on uncertainty quantification in physics-informed models via conformal prediction. It uses a neural network based surrogate (specifically FNO) as the base model and provides uncertainty via marginal and joint CP. The main innovation when compared to previous methods (Gopakumar et al 2024a) is the new score function that does not require simulator data ($Y_i$) and is instead using the surrogate itself to determine the error in the PDE. The method is tested on standard problems as well as a more advanced plasma modelling example following (Gopakumar et al 2024a).", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "# Evaluation times for CP-AER and CP-PRE:\nThank you for highlighting the need for clarity regarding computational costs. We've updated the table to separately report calibration times for both methods:\n\n| PDE  | UQ                           | L2 \\(ID\\)                    | Coverage \\(ID\\)     | L2 \\(OOD\\)                   | Coverage\\(OOD\\)     | Train time \\(hr\\) | Eval\\. time \\(s\\) | Cal\\. time \\(s\\) |\n|:----:|:----------------------------:|:----------------------------:|:-------------------:|:----------------------------:|:-------------------:|:-----------------:|:-----------------:|:----------------:|\n| Wave | CP\\-AER                      | 1\\.76e\\-05 $\\\\pm$ 4\\.40e\\-07 | 95\\.70 $\\\\pm$ 0\\.21 | 2\\.46e\\-03 $\\\\pm$ 1\\.41e\\-05 | 95\\.59 $\\\\pm$ 0\\.14 | 0:38              | 22                | 2000             |\n|      | **CP\\-PRE \\(Ours\\)** | 1\\.78e\\-05 $\\\\pm$ 4\\.61e\\-07 | 95\\.52 $\\\\pm$ 0\\.21 | 2\\.46e\\-03 $\\\\pm$ 1\\.25e\\-05 | 95\\.39 $\\\\pm$ 0\\.12 | 0:38              | 22                | **10**   |\n| NS   | CP\\-AER                      | 1\\.05e\\-04 $\\\\pm$ 6\\.58e\\-06 | 95\\.56 $\\\\pm$ 0\\.40 | 3\\.66e\\-03 $\\\\pm$ 2\\.81e\\-05 | 95\\.54 $\\\\pm$ 0\\.15 | 3:22              | 25                | 20000            |\n|      | **CP\\-PRE \\(Ours\\)** | 1\\.07e\\-04 $\\\\pm$ 5\\.18e\\-06 | 95\\.44 $\\\\pm$ 0\\.22 | 3\\.70e\\-03 $\\\\pm$ 4\\.23e\\-05 | 95\\.57 $\\\\pm$ 0\\.14 | 3:22              | 25                | **100**  |\n| MHD  | CP\\-AER                      | 2\\.20e\\-03 $\\\\pm$ 4\\.38e\\-05 | 95\\.61 $\\\\pm$ 0\\.26 | 4\\.69e\\-02 $\\\\pm$ 8\\.18e\\-04 | 95\\.60 $\\\\pm$ 0\\.27 | 5:00              | 40                | 30000            |\n|      | **CP\\-PRE \\(Ours\\)** | 2\\.20e\\-03 $\\\\pm$ 4\\.96e\\-03 | 95\\.54 $\\\\pm$ 0\\.18 | 4\\.71e\\-02 $\\\\pm$ 1\\.06e\\-03 | 95\\.67 $\\\\pm$ 0\\.22 | 5:00              | 40                | **400**  |\n\n\nEvaluation (Eval.) time refers to the time taken to evaluate the FNO. Calibration (Cal.) times refer to the time for data generation (CP-AER) of 1000 simulations or residual estimation (CP-PRE) over 1000 predictions. All evaluation and calibration times are done on a standard laptop, while the training is done a single A100 GPU. \n\nCP-AER requires substantial computational resources for calibration data generation, often through complex finite volume/element simulations that demand domain-specific expertise. CP-PRE leverages finite difference stencils as convolutional kernels, enabling GPU-parallelized computation through ML libraries with simultaneous space-time evaluation and cross-domain transferability. Our additive operator structure balances computational efficiency with statistical sufficiency (Appendix D).\n\n***\n***\n# Guidance on when CP-AER is preferred over CP-PRE:\nThank you for this important question. Each method offers distinct advantages depending on the application context:\n\nCP-AER is preferred when:\n\n•\tPhysics knowledge is incomplete or uncertain\n\n•\tBounds on physical vector fields are specifically needed\n\n•\tSufficient computational budget exists for calibration data\n\n•\tThe system has unknown or complex physics, unable to be expressed in a residual form\n\nCP-PRE is advantageous when:\n\n•\tCalibration data is prohibitively expensive\n\n•\tComplete physics knowledge is available to get bounds on conservative variables\n\n•\tComputational efficiency is critical\n\nThe primary limitation of CP-AER is its dependence on calibration data, which becomes expensive as simulation complexity grows. Meanwhile, CP-PRE only requires sufficient physics knowledge to formulate the equality constraint needed for the residual calculation.\n\n\n***\n***\n# Broader discussion of the results and their implications:\nWe thank the reviewer for this suggestion and upon much deliberation we agree that it’s more informative to provide the tables in demonstrating coverage as opposed to the figures of the residuals and the bounds. The tables (2-5) help illustrate our idea both qualitatively and quantitatively without needing any domain knowledge. The figures showing the bounds over the physical fields are indicative of the conservation laws associated with the problem and might not provide clarity to the uninitiated reader. We have the meaning of the figures in detail for the same query raised by reviewer MQ8g (under utility of method). We have restructured the paper to abide by your suggestion. Thank you for your valuable advice. \n\n***\n***\n# Primary area as Applications->Chemistry, Physics, and Earth Sciences:\nOur primary area selection reflects the paper's core contribution to physics and computational science applications. While using established benchmarks, our framework addresses a significant gap in scientific computing by providing statistically guaranteed, physics-informed uncertainty quantification for neural PDE surrogates.\nSections 5.4-5.5 demonstrate applications to fusion plasma modelling and tokamak magnetic equilibrium, showing how our method enhances reliability in high-consequence scientific domains with rigorous uncertainty bounds.", "claims_and_evidence": "The claims (as listed in lines 27-29, column 2) are supported.", "methods_and_evaluation_criteria": "The proposed method and evaluation criteria make sense for the given problem though I have some concerns over the methodology itself (see questions below). Some benchmark examples and comparisons are provided in the supplement.", "theoretical_claims": "I did not check every step in the theoretical results (Supplement A) but I believe the results follows from standard CP results as long as the residuals are exchangeable.", "experimental_designs_or_analyses": "The design of experiments closely follows previous studies, e.g. (Gopakumar et al 2024a) and similar literature.", "relation_to_broader_scientific_literature": "Yes, the key contributions are related to the broader literature. Table 2 in Supplement C provides a high level overview of how this paper fits within the broader literature and some discussion on relevant literature in both Physics-informed ML and CP is given throughout the paper.", "essential_references_not_discussed": "Not that I can tell.", "other_strengths_and_weaknesses": "Could you discuss the results in Tables 3-5 in more detail? What is the Evaluation time, particularly in the context of the two CP based methods. As far as I understand, the CP-AER method will need to evaluation the FNO and run the simulator once for each CP sample while CP-PRE does not require the simulator but requires the various differential operators to be applied to FNO. I assume the evaluation of FNO is very fast and I can see how for the examples in your paper the simulator evaluations may also be very fast as long as an optimised solver appropriate for the given problem is used. However, I do not have a good intuition of how expensive are the evaluations of the differential operators when applied to FNO but I would guess they are not much better than the simulator since they are performing ultimately the same calculations and the simulator may be highly optimised (e.g. FEM using high performance linear algebra libraries). Can you comment more on this and hence the resulting evaluation times? Can you then provide some guidance on when CP-AER is preferred over CP-PRE?\n\nThe results shown in the main part of the paper are hard to interpret. For example, looking at Fig. 4 or Fig. 8, it is not obvious whether the method does well or poorly, or how the CP uncertainty estimates could be interpreted. This is further obscured by the fact that the predictions of the model are not in the original parameter/observation space but it in the residual space, as also pointed out by the authors. A broader discussion of the results, their interpretation and the implications would make the Experimental section stronger. I would also suggest replacing some of the figures in the main part of the paper with the model comparison results that are currently given in Supplement C, e.g. Table 2 and Tables 3-5. I believe these offer more insight than the visual representations of the solutions.", "other_comments_or_suggestions": "The paper is overall well-written. \n\nLine 368: Feel - fail", "questions_for_authors": "The primary area for this paper is  Applications->Chemistry, Physics, and Earth Sciences. However, the problems shown in the paper are closely aligned with examples used in existing literature and do not solve a new problem (as far as I can tell). Is there a particular reason you chose this primary area?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper introduces a physics-informed conformal prediction framework for neural PDE solvers, providing model-agnostic uncertainty quantification without labeled data. Reviewers acknowledged the method's strong theoretical foundation and empirical validation across multiple PDE systems. Concerns about implementation clarity and computational tradeoffs were noted, particularly regarding practical interpretation of results. These issues do not outweigh the work's technical contributions, but should be addressed in revision. The framework advances reliable deployment of neural PDE solvers, warranting acceptance if there is room.
Decision: Accept (poster)