ComPhy: Composing Physical Models with end-to-end Alignment
Authors: Alessandro Trenta
Andrea Cossu
Davide Bacciu
Venue: ICLR 2026 Poster (2026)
Keywords: Learning physics
Physical systems
Partial differential equations
Systems of PDEs
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We introduce ComPhy, a multi-module approach to learn systems of PDEs by assigning one equation to each module. An alignment mechanism ensures the networks share information to solve the system together.

Abstract: Real-world phenomena typically involve multiple, interwoven dynamics that can be elegantly captured by systems of Partial Differential Equations (PDEs).   However, accurately solving such systems remains a challenge. In this paper, we introduce ComPhy (CP), a novel modular framework designed to leverage the inherent physical structure of the problem to solve systems of PDEs. CP assigns each PDE to a dedicated learning module,  each capable of incorporating state-of-the-art methodologies such as Physics-Informed Neural Networks or Neural Conservation Laws. 
  Crucially, CP introduces an end-to-end alignment mechanism, explicitly designed around the physical interplay of shared variables, enabling knowledge transfer between modules, and promoting solutions that are the result of the collective effort of all modules. 
  CP is the first approach specifically designed to tackle systems of PDEs, and our results show that it outperforms state-of-the-art approaches where a single model is trained on all PDEs at once.

Published as a conference paper at ICLR 2026
COMPHY: COMPOSING PHYSICAL MODELS WITH END-TO-ENT ALIGNMENT
Alessandro Trenta; Andrea Cossu, Davide Bacciu
Department of Computer Science
University of Pisa, Italy a alessandro. Arrent@@phduumip...is unipi.it, andrea. cossugunipi. davide bacciuuunipi.
ABSTRACT
Real-world phenomena typically involve multiple, interwoven dynamics that can be elegantly captured by systems of Partial Differential Equations (PDEs). However, accurately solving such systems remains a challenge. this paper, we introduce ComPhy (CP), a novel modular framework designed leverage the inherent physi- cal structure of the problem to solve systems of PDEs. CP assigns each PDE to a dedicated learning module, each capable of incorporating state-of-the-an method- ologies such as Physics Physics-Informed Neural Networks or Neural Conservation Laws. Crucially, CP introduces an end-to-end alignment mechanism, explicitly designed around the physical interplay of shared variables, enabling knowledge transfer between modules, and promoting solutions that are the result of the collective effort of all modules. CP is the first approach specifically designed to tackle systems of PDEs, and our results show that it outperforms state-of-the- approaches where a single model is trained on all PDEs once.
INTRODUCTION
Machine Learning is rapidly developing as a new tool to solve dynamical and physical systems, particularly for Partial Differential Equations (PDEs). These methodologies have applications in real-world problems, including weather prediction (Pathak et al., 2022), fluid modeling (Zhang et al. 2024), quantum mechanics (Mo al. (022). and molecular dynamics (Behler & Parrinello 2007 While standard supervised methods require classical solver or real sensors to gather data and train the model, unsupervised approaches can solve systems of PDEs solely from their definition. Physics-Informed Neural Networks (PINNs, Raissi et al. (2019) are considered one of the most promising approaches for physical and scientific applications. Neural Conservation Laws (NCL, (Richter- Powell et al. ) represent an alternative approach that guarantees solution with zero divergence with respect to the inputs. We formally introduce PINN and NCL. Section 2.2
When considering systems of PDEs, current approaches add a loss term for each PDE. However, this approach can be unstable and produce unwanted results. For example, PINNs often fail to converge (Krishappiyyan et al.] [222 1. Wang et al. (022) or to capture the turbulent effects that characterize fluids (Sun al., [0200] Xiang et al. 2022). Some workarounds require reweighting the loss at each epoch (Wang ||22 Xiang al. 2022 McClenny & Braga- 2023), enforcing boundary conditions Sun et al. 2020). or dividing the spatial and temporal domains into subdomains (Jagtap al. 2020 Kharazmi et al. 2021 Meng et al. 2020 Krishnapriyan et al. 2021). However, most of these approaches are not specifically ly designed for systems of PDEs, or only focus on a specific one.
In this work, introduce ComPhy (CP), multi-module approach to solve systems of PDEs by allocating one learning module for each PDE and ii) leveraging the physical structure of the problem align modules that learn the same physical variable(s). CP be applied to any system of PDEs, and the learning modules can be chosen among any of the existing approaches (we experimented with PINNs and NCL). Intuitively, each CP module solves simpler problem than the overall system, while our end-to-end alignment process integrates the knowledge from all modules to devise a collective
"Corresponding author.

---

Published as a conference paper at ICLR 2026
AD
PINN
iu
Op
+p Syy
APi By
8
[uu u vuu Vp
0(0) u(O)
ALIGNMENT
NCL PINN V u =
Q
a(0) u(0)
Figure 1: Visual representation of ComPhy for the davier-Stokes equations. Each module (in this case, PINN or NCL) learns one PDE. The derivatives of each module, computed with Automatic Differentiation (AD) (Baydin al. 2018), are used to optimize each PDE and to align the other module. CP only uses the first module for inference as it predicts all the variables of the system.
solution. Practically, the alignment encourages the outputs of different modules representing the same physical quantity to be similar. We propose several alignment losses, highlighting the role played by the derivatives in transferring physical knowledge between models (Trenta et al. 2026).
While the training phase involves all learning modules, during inference CP uses the minimum number of modules needed to predict relevant variables. Most of the time, one module sufficient. We show that CP achieves fate-of-theearr performance, outperforming alternatives where a single learning module is trained on all the equations of the system.
methodology to solve systems of PDEs by allocating one module for each equation while aligning modules that learn the same physical information. (ii) We introduce and evaluate three different alignment approaches. (iii) We evaluate several CP configurations on case studies involving systems two and three PDEs. We achieve state-of-the-art results in all cases. We analyze the backpropagated gradients through the modules, showing why CP trains better than PINNs. (iv) We evaluate CP on two challenging real-world physical systems with and equations, achieving ofte-of-theear results.
Our contributions are: (i) ComPhy,
COMPHY
We introduce ComPhy, our end-to-end and modular approach for solving systems of PDEs. Our approach involves: (i) learning each PDE with single module while (ii) constraining and refining the learning process of each module through shared knowledge and exchange of physical information. In particular, the second point is implemented via alignment loss between the modules.
2.1 METHODOLOGY
To set the stage, let us consider points æ sampled from a spatial domain 52 C R. with time in the interval [0, T]. We call 0Z the boundary of the domain Ω. We define system of N PDEs by
F [u] = 0, Bluj b(t,
(PDE i) for -= 1,,..,,
(BC)
(1)
(0, (
(IC)
where utt, .)) is the solution to the system of PDEs. A necessary condition for the existence
and uniqueness of solution (Evans 2022) is to provide the initial condition (IC) u(0, æ) g(x), as well as the boundary condition (BC) B u = b(t, x). The PDEs must be satisfied for every (t, E [0, Ω the IC for every Ω while the BC for every (t, x)) [o, T| x 02.
State-of-the -h-thhrrrt models like PINNs may suffer from optimization problems when multiple PDEs involved (Sun al. 2020 Wang et al. 2022). due to the ill ill-conditioning and the complex loss landscape of PDE residual terms (Krishnapriyan et al. (02I). which can lead to strong imbalances in the gradient norms of the different loss terms (Wang al. 2022). ComPhy avoids this issue
2

---

Published as a conference paper at ICLR 2026
by using different modules to optimize the different PDEs of the system separately. This makes each optimization problem simpler than solving the whole system, while cross-PDE information is acquired through alignment. In the CP modular architecture, each network takes as input the same coordinates (t, x) and predicts the relevant variables in its equation. Modules can sometimes solve more than one equation, but no module solves the whole system. Each module learns the IC, BC, and the PDE(s) of interest through the weighted sum of the respective loss terms:
m -BBCCBC + Aichic + APDE SyDEE
=)ncllu(t, x) b(t, x)) [[A((0.7)x00)) + Arcll x) - g ((=)) + Appn:
The loss term LYDEi depends on the particular choice of the learning module (Section 2.2 shows the CP architecture for the Navier-Stokes PDEs, which we discuss Section 2.3
Figure
Alignment. If the PDEs in the system not all satisfied at once, the solution to the system might not be unique (Evans) 2022). Hence, we need to provide each module with information about the state of the others, least on the parts with the same physical information. This is the role of our alignment process.
Let be the subset of shared variables predicted by two different modules. The alignment loss fosters consistency between the predicted v E v by both models. We focus on 2 alignment losses acting on the prediction of made by module (6,), and by module j (ěz):
SOB: Calign 1j = 0 - - + [Jeē; - Jú, 12:
DERL: La = Jo; - J
where Jé; is the Jacobian of the predictions for a of module with respect to all its inputs. SOB aligns the Sobolev norm 2 (Maz ya 2011] while DERL aligns the 2 norm the derivatives only. Incorporating derivative information improves the convergence (Czarnecki et al. 2017) and imposes stricter resemblance between the models' outputs. In physical systems, derivatives describe the system's 's dynamics and provide sufficient information to determine its evolution. Furthermore, they are the most effective target for transferring and distilling physical information between models TTenta et al. 2026). To substantiate the claim, we also consider the OUTL alignment loss 0/0; 6 | as well, which is based only on the L distance of the outputs. Given our discussion, we expect this perform much worse than the others. See Appendix [B.2 (B.2 for more details on the alignment losses.
The complete loss of the CP model can be expressed as
N Σ 1=1
CCp
Aalign Calign ij +
Lmoduke i:
where the sum on j spans each pair of modules sharing at least one variable. All L losses are implemented using the empirical version, that is, the mean squared error on collocation points. In reality, we see that it is sufficient to consider a simplified version of this loss, which includes only the alignment terms between the inference module and the others. For more details, see Appendix C. Inference. CP trains all the modules end-to-end. During inference, CP uses only one or a subset of the modules. In particular, is sufficient to include the minimum number of modules required to predict all the variables present in the solution. This often saves a considerable amount of inference time, as even a single module is often sufficient to predict the output.
2.2 LEARNING MODULES
We now introduce the learning modules we adopted for our experiments with CP. We consider existing the models that are widely used to learn physical systems. As we demonstrate in our experimental sections, Sections 3 and 4 optimizing all equations within a single module is suboptimal for learning challenging systems of PDFs due to the different scales and competing objectives of the losses. Our model, based on aligning specialized modules, overcomes these issues and outperforms all baselines.
Physics Informed Neural Networks (PINNs). PINNs Raissi et al. 2019) are among the best- performing and most used models for learning solutions to PDEs. They be constructed upon
3

---

Published as a conference paper at ICLR 2026
any differentiable architecture. As commonly done in the literature (Raissi al. 2019 Lau et al. 2024). we consider multi-layer perceptrons (MLPs). PINNs take as input the coordinates of a point (t, Φ) and output the prediction of the solution u(t, x) to the PDE. To solve the problem, PINNs try to impose the PDE soft target for the MLP itself. Partial derivatives of the MLP du 0ú 8t are calculated via Automatic Differentiation (AD. Baydin et al. (2018). Then, the residual on the PDE is calculated using these derivatives Flü]. Finally, the L loss is used minimize this residual, together with the IC and BC, obtaining the full loss:
PINN || m) b(t, 1) + || x) - g(=) |F
(5)
Neural Conservation Laws (NCL). Being able to output conservative field by design is the key feature of NCL (Richter-Powell et al. (2022)). To achieve this, the MLP does not directly output the field but parametrizes antisymmetric matrix A. Then, the divergence operator @f
is applied row-wise, obtaining the final vector u with items u, div.A..). This parameterization ensures that V uu- = div(u) = 0, so the field is divergence-free by design. Richter-Powar et al. (2022) consider all input and output variables in this process, calculating the divergence with respect to the entire input vector (t. In our case, we generalized their approach to produce divergence-free fields with respect to any subset of the input and output variables, as long as they are of the same size. See Appendix B. for more details and examples.
div(f) - 2i=1 Sæ;
2.3 A PRACTICAL EXAMPLE: NAVIER-STOKER EQUATIONS
Having introduced all the components of CP, we now go through a practical example of its application. We consider the 2D Ntvier-Stokes equations:
Hu ΦE
[
Au Vp
momentum eqn.
(6)
u = 0.
incompressibility eqn.
To solve this system with CP, we consider two modules. The first one PINN and learns the momentum equation in 6 predicting (û1, P1). Its specific loss components in equation 2 are
LBc || u [| -p1
Lic
|| (0, x)) u(0 1 || (0, x) p(0, 2)
Oû;
LPDE
+ ú, i Du
V
Øt
The second module, which is either PINN or an NCL, is related to the incompressibility equation in and predicts only u. If the module is PINN, the loss components are
LBC = Ú0 2
Lic2 = Ú x) - u00, 11.
LPDE = I/V u
(8)
while 2 is omitted this module instead an NCL, as is satisfied by design. To define the alignment loss, we identify the common variables of the two modules, that is Hence, this loss is defined as Lalign رش 12 During evaluation, we need to predict all variables that define the problem: u and p. Since the first module predicts both, it is the only module required to evaluate CP.
2.4 THEORETICAL CONSIDERATIONS
Theoretical results on the convergence of physical models (like PINN and NCL) to the true solution remain open problem. Some insights exist only on a limited class of single PDEs (Yeonjong Shin et al. 2020). Systems of PDEs represent a harder problem, as PINNs often fail to learn the underlying solution. This is especially true for the Navier Stokes equation, where the turbulent behavior of fluids represents an issue for these models (Sun et al., 2020 Xiang et al.) 2022).
Since the Sobolev distance is a widely used tool for the convergence, existence, or uniqueness of solutions to PDEs (Evans 2022). we expect our derivative-based alignment (SOB and DERL) to outperform the OUTL alignment.
Intuitively, if the SOB/DERL alignment loss always zero, the aligned modules would predict the same solution. This would make our modular CP equivalent to a single learning model that has access to the overall system. However, in this case CP would enjoy simpler optimization problem,
4

---

Published as a conference paper at ICLR 2026
Table I: Summary of the tasks we consider in the experiments. The last column indicates the different CP configurations, in terms of the modules, tested in the specific experiment.
Experiment
Equation
Description
CP Modules
Navier-Stokes
(TG.M)
[D a du Up
Time dependent
PINNS (2<PINN)
Thyfor-Green vortex
PDEs
PINNINCI.
(KT.)) [D (T..V
Navier-Stokes
a du
Up-U
Time independent
2 (2sPINN)
Kovasznay flow
(KEDV
u 0
PINN+NCL
dp
de
(A.P)
0r
+
Bty
=0
Time dependent
PINNs (3xPINN
Acousties equations
0u dp
du
1
PDBs
NCLs
(A.V%)
(A.
p.y
Navier Stokes Euler (NS-Enler) Gas equations
(EG.C) Up
+ (pw) 0 (EG.I)
V
Time dependent
PINNS (2<PINN) PINN+NCL
+ [D_uju Ve =0
PDES
(2NNCI)
(BG M)
3 PINNS (3PPINN)
(MHD.C)
(pu)
(MHDD.) V
OV
PINNS (2<PINN)
MaguetolTldrroDyymmiic (MHID)
(MHD.M)
[Dauju
= -V B) x B Vp
Time dependent PDEs
2xPINN
PINNS (3PPNNN)
BB
PINNS (1APINN)
(MHD..)
(MHD.S) de
each module would only focus equation. Furthermore, with the simplified alignment loss, each module optimizes the IC, BC, one PDE loss, and one term of the alignment loss on average, for total of 4. Instead, model optimized Hence, from PDEs onwards, CP has this further advantage. In Section 3.4 we analyze how the different loss terms backpropagate through the different layers of the modules, showing that those of CP have much more similar scales than those of PINNs, improving the quality of training steps. In Appendix G we provide further intuitions the differences between alignment and PDE residual losses, as well as a more explicit example on the transfer of physical constraints between modules.
the full system of PDEs involves
loss terms.
3 EXPERIMENTS, CASE STUDIES
We empirically validate the effectiveness of CP on several physical systems of PDEs. We work with systems of 2, 3 and 5 PDEs. Parvvi[] summary of the tasks and the related CP configuration. In this Section, we focus on the Navier-Stokes and Acoustic equations as the first case studies. At the end, we also provide an empirical explanation of why CP outperforms state-of-the-ar approaches by analyzing the gradients backpropagated through CP modules. Section ! evaluates CP on two challenging real-world systems: NS-Euler Gas and MagnetoHyydroyynamic (MHD) equation
We compare CP with four main approaches: PINN (Raissi et al., 2019). PINN with gradient based reweighting (PINN+Grad, Wang et al. (2023), ) PINN with RAR point resampling (PINN+RAR, Wu et al. (2023)) and NCL (Richter-Powe Powell et al. 2022). Further details are available in Appendix These models are implemented as single model that accesses the information about all PDEs. This is done by replacing the last term Equation 5 with Σ [|FrDs: [u]]l We remark that NCL satisfies divergence-free ffre equations by design, without requiring any loss term. We remark that NCL a generalization of the original approach of Raissi et al.] (2019) for divergence- freee fields. For completeness, also test CP with gradient-based reweighting.
Following Wang et al. (023), we train all models for 100, 000 steps with batches of 1000 points randomly sampled in the domain. We the Adam optimizer (Kingmal 2014) with initial learning rate of 0.001, decaying every 2000 steps by factor of 0.9. Model architectures are in Appendix CC3
To measure the accuracy of the predictions, evaluate all approaches with two metries: (i) |/u ûlI2 (L²-err), which measures the 2 distance between the true solution u and the estimate á. (ii) (,, l ui (max-err), which measures the maximum absolute error in the domain.
'The code reproduce the experiments available http://ghhhbbccmm//////httyyCym
5

---

Published as a conference paper at ICLR 2026
Table Results for the Taylor-Green vortex, Kovasznay flow, and Acoustics experiments. L errors of the model predictions with respect to the ground truth, and the maximum ertors in the domain.
Taylor-Green Vortex
Kovasznay Flow
Acoustics
Model
Aligument
L7-en:
L--err.
Model
Aligument
L'-err. x10-3
x10-5
x10-3
x10-6
x10-3
x10-
PINN
3.677
9.374
10.83
5.203
PINN+RAR PINN+Grad
5.235
6.160
11.48
8.062
PINN-RAR
120.6
20.05
4.245
9.719
7.471
4.264
GINN+Grad
10.05
3.391
NCL
2,830
5.947
3.014
1.398
NCL
5.243
1.552
CP-2xPINN CP-2XPINN
OUTL
6.028
9.945
20.47
4.972
CP-3sPINN CP-3XPINN
OUTL
11,43
8.437
SOB
5.222
8.027
8.641
6.195
SOB
7.611
2.148
+Grad
4.687
6.250
7.695
4.207
+Grad
5.726
1.685
CP-2XPINN
DERL
3.619
6.378
7.221
4,996
CP-3XPINN
DERL
7.449
2.171
+Grad
+Grad
5.718
1.737
CP-PINN+NCL CP-PINN+NCL
OUTL
6.024
11,55
13,30
4.515
CP-3NNCL CP-3xNCL
OUT.
6.008
4.549
SOB
4.764
9.215
6.867
3.619
SOB
5.172
1.521
+Grad
4.386
8.418
3.663
1.277
+Grad
2.727
1,090
CP-PINNANCL
DERL
3.247
7.996
5.422
3.468
CP-SANCL
DERL
5.165
1.547
+Grad
+Grad
PINN
PINN-Grad
NCL
CP-3XPNN (DERL+Grad)
CP CP-3KNCL (DERL+Grad)
0.15
00
0.55
Figure 2: Acoustics equation: errors for PINN, PINN+Grad, NCL, and the best CP models (one for architecture) at = 0.16. The lowest errors are in blue.
3.1 TAYLOR- VORTEX
Our first experiment involves the 2D time-dependent incompressible Navier-Stokes equations, namely the momentum (TG.M) and incompressibility (TG.I) Equations Table [ We work in the domain [0, T] x [0 with time f E [0. 10) and a viscosity coefficient of v = 0.1. The analytic solution of the Taylor-Green vortex (Chorin 1968), used calculate the reference solution, is given Appendix
D.1 The inference module for CP is the one learning equation (TG.M), as it predicts all the variables.
Table shows the numerical results for this experiment, while Appendix [.J]shows the prediction error for different values of t. We see that CP models aligned with DERL outperform all variations of PINNs. CP-PINN+NCL with DERL alignment and Grad reweighting performs best the 1 metric, even when compared to NCL. This is remarkable, since NCL. has the built-in advantage of bypassing the divergence-free equation, while CP has to transfer it across modules. Additional results available in Appendix D.I
3.2 KOVASZNAY FLOW
We consider the 2D lime-independent Navier-Stokes equations made of the momentum (KF.M) and
incompressibility (KF.I) equations in T We work in the domain [-1, 1] 0-0.5. 1.5 with
viscosity coefficient of v = 50
1 The analytic solution of the Kovasznay flow Drazin & Riley 2009
which represents the flow behind two-dimensional grid, is given in Appendix D.2 The evaluation
model for CP is always the one for equation (KF.M), as it predicts all the variables. Numerical results
are provided in Table The complete of results and figures is available in Appendix D.2 None of the PINN variants matches the performance of the CP models aligned with either DERL or SOB. Interestingly, OUTL alignment delivers worse results, supporting our hypothesis that derivative-based alignments are more effective. Due to its implicit advantage of modeling divergence- freeeeeeeeeeeeeeeeeeeeeeeee fields, NCL shows a competitive performance. However, CP with gradient- based reweighting obtains the best performance possible, between 2 and times better than PINN variants. We conclude that within these Navier-Stoker case studies, ComPhy alleviates PINN optimisation issues and performs comparably to, if not better than, NCL. However, CP is far more general than NCL,, and can be

---

Published as a conference paper at ICLR 2026
Layer
Layer 2
Layer
Layer
PrQquetty meeqengy Requero
104
201
Grad Viiue
vely
algnment
pres
(a) CP-3xPINN (DERL)
Layer
Layer 2
Layer
Layer
Rrquency
11
02
0.00
0.02
-0.02
0.00
0.02
0.00
0.02
-0.00
Gradient Value
vely
Vel
(b) PINN
Figure 3: Acoustics experiment gradient histograms. Each plot contains the histograms for the distribution of the gradients propagated at each layer the CP-3xPINN and PINN models at the beginning of training, similarly to Wang et al. (2021). Different colors are for different losses.
directly applied to systems without divergence-free free components. As OUTL consistently showed worse results, for the remaining experiments, we only work with derivative-based alignment.
3.3 ACOUSTIC EQUATIONS
We now consider set of three equations, namely equations (A.P), (A. and (y. in Table T solved using CP model with three aligned modules. As divergence-Fe form (see Appendix [B.1), respectively for t, x, for the first equation, t, x and t. y for the second and third. Hence, we can employ (A.P) contains all relevant variables, the first module is the one used for evaluation.
see, all three equations are in the
CP model with three NCL modules. Since equation
We work in the spacetime domain of (t, x. y) e [0, 0.24] x [-1, 1]3, with density م and bulk coefficient K = I. Appendix D.3 provides the initial conditions, while boundary conditions are null everywhere. The reference solution calculated using the Clawpack software package (Clawpack)
Development Team 2024 Mandli et al. 2016). Further details can be found in Appendix D.3
We report the numerical results in Table while Figure 2 shows the prediction error in the domain at t= 0.16. For additional Figures, see Appendix D.3 Even when using 3 different modules, CP is clearly the best approach for both our metries. In particular, CP with 3 NCL modules performs best. None of the PINN variants fall close to our best model, not even when adding resampling or adaptive weights. CP, instead, benefits from gradient-based reweighting, showing a 2x/3x improvement with respect PINNs and NCL. This experiment suggests that CP and derivative-base alignment scales well to a larger number of modules. Section scales CP to even larger systems.
3.4 GRADIENT ANALYSIS
To get a better understanding of why CP outperforms PINNs, we conduct an empirical analysis on the distribution of the backpropagated gradients different layers of the CP modules. Given w, the set of weights in layer j. for each loss term LBC, Lic, LPDE; of equation 2 we collect the gradients dc and plot the corresponding histogram grouped by layer (Figure
A similar analysis for
Buyy
PINNs is discussed in Wang et al. (202I) The authors show that when PINNS are optimized correctly, the gradients of the different losses tend to be evenly distributed within each layer. However, if the distributions differ significantly, implies that one loss function is dominating, which negatively affects the model's performance.
We replicate this study for CP on the Acoustics experiment of Section 3.3 which involves 3 equations and 3 modules. Figures 3a and 3b show the distributions of propagated gradients for CP 3xPINN
7

---

Published as a conference paper at ICLR 2026
Table Numerical results for the Euler gas and MHD equation experiments. We report L emrors of the model predictions with respect to the ground truth, and the maximum errors in the domain.
(a) Euler gas equation experiment.
(b) MugnecooHyydooDynamiis experiment.
Model
Alignment
L2-emr x x10-3
max err
Model
Alignment
L2-emr. x10
max em
x100
x10
PINN
1.712
5.364
PINN
_
1.967
12.74
PINN+Grad
3.319
7.127
PINN+Grad
1.657
9.164
NCL
1.690
4.912
NCL
1.975
12.59
CP-2xPINN
SOB
1.296
3,082
CP-2xPINN
SOB
1.599
9.324
CP-2xPINN DERL.
1.401
3.149
CP-2xPINN
DERL.
1.608
9.280
CP-PINN+NCL
SOB
1.487
2.721
CP3xXPNN
SOB
1.535
8.528
CP-PINN+NCL
DERL
1.468
2.673
CP-3xPINN
DERL
1.535
CP-2xNCL
SOB
2,029
4.311
CP-2XPINN+NCI
SOB
1.601
8.913
CP-2xNCL DERL
1.396
4.183
CP-2XPINN+NCI.
DERI.
1.607
8.934
CP-3xPINN
SOB
1.382
3.021
CP-4XPINN
SOB
1.560
8.764
CP-3xPINN
DERL
1,462
2.899
CP-4xPINN
DERL.
1.567
8.771
30
PMN
P
NCL
CRAINNANC IDERL
CP 2PPNPINN (DERU)
(DERL)
CP IDERL
CP CP33PPNN (BERL).
28
20
5 10:00
D0P
LS5 010
L:-00
0,
L0:00
10
Figure 4: NS-Euler gas equations experiment: prediction error in the domain at / 0.3 for the PINN, NCL, and the CP model (best for each architecture).
and PINN, respectively. The plots clearly highlight that combining three modules makes the gradient distribution more aligned across layers. As we have seen in Section 3.3 this also entails better performance. We provide the same analysis on the Kovasznay flow experiment in Appendix D.2 Both studies support the fact that the modular approach ComPhy results in an casier optimization problem when compared single PINN trained on the entire system of PDEs.
EXPERIMENTS: REAL-WORLD SYSTEMS
4.1 NS-EULER EQUATIONS
We now challenge CP with two real-world systems: the NS-Euler and the MagnetoHlydroDynamiss equations. The NS-Fuler equations for incompressible gas are given by the conservation (EG.C),
incompressibility (EG.I), and momentum (EG.M) equations in Table T] where is the density, u the velocity and p the pressure. Our setup is the same as the in Richter-Powell
(2022): :
consider the 2D torus and work with the spacetime domain (t, v) 0, 3 [0, 1 Appendix D.4
provides the initial conditions and further details. The reference solution available and collected from Michter-Powell et al. (2022), and was calculated via the Finite Element Method (Anderson al. 2010)All models are trained for 600, 000 steps with batches of 1000 random points at each step. The optimal hyperparameters and architectures are taken from Kichter-Powwll et al. (2022). On the NS-Euler equations we test CP configurations: CP-2xPINN, CP-PINN+NCI. CP-2xNCI, and CP-3xPINN, For the first three, the first module (cither PINN NCL) learns both equations (EG.C) and (EG.M), while the second learns only (EG.I). Since all relevant variables are in equation (EG.M), only the first module is used for evaluation.
Numerical results are provided Table 3al while Figure shows the prediction error at = 0.3. Appendix D.4 shows the model errors for different values of t [o. 21. While Richter-Powad et al. (2022) clam that the PINN fails to converge to the true solution, we did not encounter the same issue. Our PINN solve the system some extent. CP achieves state-of-the-ar results when using both
8

---

Published as a conference paper at ICLR 2026
CH-SenNAN BUeno
NiL IteHi
Figure 5: MHD experiment: prediction error in the domain = 0.5 for the PINN, NCL. and the CP model (best for each architecture).
DERL and SOB, with the latter struggling only when 2 NCL modules are employed. This highlights the role played by the outputs' derivatives in the alignment: a loss computed only from outputs' derivatives (DERL) is much more effective than a loss computed also from the output (SOB). We conclude that CP learn challenging real-world system of PDEs effectively, using different combinations of NCL and PINN models.
4.2 MAGNETOHYDRODYNNANN
Finally, we scale our model to system of 5 equations, the MHD equations in Table which represent the motion of plasma with combination of the compressible Navier-Stokes equations for fluids with the Maxwell equations of Magnetism (Gruber & Rappaz 1985). These are composed of the continuity (MHD.C), momentum (MHD.M) (M), state (MHD.S), induction (MHD.I), and Gauss (MHD.G. G) equations. We consider a similar setting to Gopakumar et al. 20255), with domain 0. 1 and time E [o, 0.5]. The reference solution is calculated with the finite volume method (Ferziger & Peric 2001). Since we are interested in pure scaling, we consider combinations of 2, 3 and PINNs for our CP models, as well as NC We compare a PINN, PINN+Grad, and an NCL. with equation (MHD.G) satisfied by design. Architectures are available in Appendix C.3 Models are trained for 100, 000 steps with the Adam optimizer, and batches of 1024 randomly sampled points.
We report the results in Table 3b with Figure SJohiiing the prediction error at = 0.5. Even with 5 equations, CP outperforms PINN and NCL by a large margin with any the tested combinations of modules, including that with 4 modules. We report additional results and Figures in Appendix D.5
4.3 DISCUSSION AND FURTHER EXPERIMENTAL ANALYSIS
On the choice of the CP modules. Apart from the MHD equations, where we tested the scaling capabilities of CP to a higher number of PDEs, we experimented with every possible combination of modules, keeping in mind that the NCL module requires divergence-free equation to be applied. Our results support the fact that, when possible, a combination of modules that includes NCL is always beneficial. Similarly, since most physical systems present 3 equations, the most common and effective choice is to use 2 or 3 modules, respectively.
On the strength of the alignment process. In Section 2.4 and Appendix we provided theoretical understandings on why the alignment mechanism allows for effective transfer of physical constraints between modules. When the alignment loss converges zero, the corresponding modules behave the same function in the Sobolev space 1,2 (Trenta et al. (2026). To empirically validate these claims, we plot the curves of the alignment loss, the prediction loss, and the PDE residuals that are not actively trained on the inference module. We also empirically measure the distance and relative discrepancy between the predictions of different modules. The results, available in Appendix [E] confirm theoretical claims and the strength of the alignment mechanism: modules predict very similar functions, while prediction and PDE losses are correlated with the alignment loss.
Ablations. To test the robustness of ComPhy in different setups and conditions, we performed ablation studies available in Appendix F (i) We tried different values for the alignment coefficient Aalign- showing that the naive choice of Aaignn = works well on average, with possible improvements with highly specialised tuning of this parameter. We remark that the Gradient-based reweighting scheme used in Section alleviates the burden of hyperparameter selection. (ii) We incremented the number of units in the MI.Ps of the PINN and NCL baselines to match the parameter budget of

---

Published as a conference paper at ICLR 2026
CP during training. This leads to fewer parameters used by CP during inference, as it employs only one module at that state. Even so, CP performs better, showing the importance of modularization for these applicatiions. (iii) We run the same experiment multiple times with different seeds to test the robustness to initializations and stochastic gradient descent. The results remain strong, as the CP models have the best mean results by far, and lower standard deviation.
RELATED WORKS
PINNs Raissi et al. 2019) are widely used to solve physical problems such as fluid dynamics (Sun et al.) 2020). However, PINNs suffer from optimization problems, such as imbalances in the gradients of the different loss terms (Wang et al., [2022 or the ill-conditioning and complex loss landscape
of PDE residual terms (Krishnapriyan et al. 2021). PINNs can also fail to reach a minimum of the loss (Sun et al. 2020). Recent approaches try to alleviate these problems with resampling strategies that explore regions with higher PDE residuals while retaining enough points in well-optimized regions (Daw et al. 2023), or by using a variational formulation of the PDE residual term in small regions around collocation points (Wu et al. 2024). Other works adopt strategies to align the gradients at different gradient descent steps using second-order corrections (Wang et al. 2025) which can be computationally demanding, or to align the autoregressive updates to the underlying true dynamics (Zhu et al. 2025) for data-driven models. Instead, we work in an unsupervised setting, and our alignment imposes on the different modules to share information on their respective physical constraints. Even though many other solutions, such as resampling or reweighting schemes, have been proposed (Sun et al.) 2020; Wang et al. 202T; Xiang et al., [022] McClenny Braga-Neto 2023 Zhao |021 lau et al. 20024), none of them is tailored to tackle systems of PDEs.
Ad-hoe architectures solve the optimization issues by satisfying physical properties by design. NCL (Richter-Powell et al. 2022) ensures that the output of the model is divergence- while Torres ctt al. (2024) proposes a divergence-free normalizing flow. However, such solutions only work for a given physical property and still require enforcing additional ones with alternative approaches, like PINN. Finally, our CP is the first attempt at solving systems of PDEs with a modular, end-to-end approach where modules are aligned to transfer physical constraints with each other.
Neural Operators (Kovachki et al. 2023) are deep models that learn mappings between functional spaces. DeepONets (Lu al.. 2001 try learn representations of the input function and combine them on the output function inference points. Graph Neural Operators (Li et af.) 2020), instead, transform functions via local convolutions with parametrized kernels, which are learned during training. Finally, Fourier Neural Operators (Li et al. (02I) use Fourier transforms to simplify and speed up the calculation of such convolutions while acting on the function globally. While these methods are powerful for supervised learning tasks and generalization to different initial conditions or parameters, they are more computationally demanding and serve a different task, as they learn mappings between functions. Furthermore, they need simulated data, and they are not suitable for the unsupervised setting of PINNs.
6 CONCLUSION
We introduced CP, a modular approach for solving systems of PDEs where the final predicted solution is the result of the collective effort of the modules. The modules communicate with each other during training via alignment, a process that encourages the predictions of the same variable made by different modules to be similar. In particular, we consider a derivative-based alignment, and we show that it provides state-of-the-art results on several systems of PDEs, including challenging world tasks. We further analyzed why CP optimizes the thysics-informed losses better than PINNs, finding that their propagated gradients have more similar scales. While NCL is competitive in the case studies, the advantage of having one equation less to be optimized fades away in real-world systems of three or more equations, where CP scales much better.
Future works could consider different optimization protocols, such as alternating module optimization steps with aligmment steps, as well as aligning higher-order derivatives to capture richer dynamics of the physical system. CP represents an initial step toward developing modular architectures for scalable learning of physical systems. Progress in this direction could lead to dynamic models capable of adapting to new initial conditions, as well as incorporating or removing modules over time to accommodate changing environments.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work has been supported by EU-EIC EMERGE (Grant No. 101070918).
ETHICS STATEMENT
The research conducted in this work completely adheres the ICLR Code of Ethics. This work does not involve human subjects, nor does it provide harmful insights, discrimination/bias concerns, or any other potential ethical concems.
REPRODUCIBILITY STATEMENT
We provide all the necessary information reproduce the results in the paper in the main text and in the Appendix. The methodology is described in Section while experimental settings, along with training steps, optimizers, and model architectures, are described in Section 3 and Section Further implementation details are available in Appendix[[] We additionally provide all the code necessary to reproduce the experiments, from data generation to model definition and training, as well as evaluation, as Supplementary Material. The code is also available at https::/github.ccmm Alex Taw.eeexhit//CoPlll
REFERENCES
Dale Anderson, John C Tannnhill, Richard H Pletcher, Ramakanth Munipalli, and Vijaya Shankar. Computational fluid mechanics and heat transfer. CRC Press, Boca Raton, FL, edition, December 2020.
Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: survey. Jourmal of Machine Learning Research, 18(153):1 43, 2018. URL URL http:://mma...cccaaaaatt////////
html
Jörg Behler and Michele Parrinello. Generalized neural-network represtttttoontnnnntnntttttooont of high-dimensional potential-energy surfaces. Phys. Rev. Lett. 98:144401, Apr 2007. doi: 10.110//hysRRewLettt... 146401. http://t////./....././//////////
Alexandre Joel Chorin. Numerical solution of the navier-stokes equations. Mathemattics Computa- tion, 1244):775576,, 1968. ISSN 00255718. 10886842. URL stable/2004577
Clawpack Development Team. Clawpack software, 2024. URL Version 1.1..
Wojciech M. Czarnecki, Simon Osindero, Max Jaderberg, Grzegorz Swirsze7, and Raz- van Pascanu. Sobolev training for neural networks. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Ad- vances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.. 2017. URL httpss://proceedinngggg//// http::ssssiiirrsreeseeeoocc flep/580611000000cc33eeee2222222
files/ ttt.e//ppppr/211....
Arka Daw, Jie Bu, Sifan Wang, Paris Perdikaris, and Anuj Karpatne. Mitigating propagation failures in physics-informed neural networks using relain-reammple-reless (r3) sampling, 2023, URL
https::/aaxxv.orr///ass
P G Drazin and N Riley. London mathematical society lecture note series: The Navier-stokes equations: A classification of flows and exact solutions series number 334. Cambridge University Press, Cambridge, England, November 2009.
Lawrence C Evans. Partial differential equations. American Mathematical Society, Providence, RI, March 2022.
Joel H Ferziger and Milovan Peric. Computational methods for fluid dynamics. Springer, Berlin, Germany, 3 edition, November 2001.
11

---

Published as a conference paper at ICLR 2026
Vignesh Gopakumar, Ander Gray, Lorenzo Zanisi, Timothy Nunn, Daniel Giles, Matt Kus- ner, Stanislas Pamela, and Mare Peter Deisenroth. Calibrated physics-informed uncertainty quantification. In Forty-second International Conference on Machine Learning, 2025. URL
htttps://cpenreviiee net www.::/wwonnoooooomooooooouuuaa
Ralf Gruber and Jacques Rappaz. The Ideal MHD Model. pp. 34-41. Springer Berlin Heidelberg, doi: 10. 10.1007/978-3-64--867777
Berlin, Heidelberg, 1985, ISBN 978-3-642-886788-88
URL
https://doi. org/10 112.0.1.00000 1001100779777 3-2-
Ameya D Jagtap, Ehsan Kharazmi, and George Em Karniadakis. Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. Comput. Methods Appl. Mech. Eng., 365(1 13028): 13028, June 2020.
Ehsan Kharazmi, Zhongqiang Zhang, and George E M Karniadakis. hp-VPINNS: Variational physics- informed neural networks with domain decomposition. Comput. Methods Appl. Mech. Eng., 374 (1 1133547):113547, February 2021.
Diederik Kingma. Adam: A method for stochastic optimization. arXiv preprint arXiv: 1412.6980 2014.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Reseurch, 24(89):1-97, 2023. URL
fttp::/jjmi= org/ papers/
wtteers/v244222 wtt///24/2222222 1524 html
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. In M. Ran- zato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 26548-26566 Curran Asso- ciates, Inc., 2021, URL
paper. files/papes 2ttp////1.c.42200011000002227720
Gregory Kang Ruey Lau, Apivich Hemachandra, See-Kiong Ng. and Bryan Kian Hsiang Low. PIN- NACLE: PINN adaptive collocation and experimental points selection. The Twelfth International Conference Learning Representations, 2024. BRLhttt:: /ttp::///op/eeeeeeee net forrmm?rmm
1d=GzNaCp6Vcc
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations, 2020. http/t:tta//axiiii www.aaawwww...ooaaaaooa.aaeooeoe abs /tw.//00..... 03488533.0443
Zongyi Li. Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL
https:://peerreviee net
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem operators. Nature Machine Intelligence, 3(3):218-22) March 2021. ISSN 2522-5839. doi: 10.1038/ $42256-021--00002-5 URL /ttp:://dxxxxo org 10.0 10301008///0 1el.1038/422255 0220-00302-22
Kyle T Mandli, Aron J Ahmadia, Marsha Berger, Donna Calhoun, David George, Yiannis Had- jimichael, David Ketcheson, Grady Lemoine, and Randall J Leqeqe Clawpack: building an open source ecosystem for solving hyperbolic pdes. PeerJ Computer Science, 2:e68, 2016. doi: 10.7717/peerj-os.66
Vladimir Maz'ya Sobolev spaces. Grundlehren der mathematischen Wissenschaften. Springer, Berlin, Germany, edition, February 2011.
Levi D McClenny and Ulisses M Braga-Neto. -dlf-adaptive physics-ifformed neural networks.
Journal of Computational Physics, 474(111722):1 February 2023.
12

---

Published as a conference paper at ICLR 2026
Xuhui Meng, Zhen Li, Dongkun Zhang, and George Em Karniadakis. PPINN: Parareal physics- informed neural network for dimeddependent PDEs. Comput. Methods Appl. Mech. Eng., 370 (113250):113250, October 2020.
Yifan Mo, Liming Ling, and Delu Zeng. Data-driven vector soliton solutions of coupled nonlinear schrödinger equation using a deep learning algorithm. Phys. Lett. A, 421((17739)) Januray 2022.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadeneshelii Pedram Hassanzadeh, Karthik Kashinath, and Animashree Anandkumar. Fourcastnet: A global data- driven high-resolution weather model using adaptive fourier neural operators, 2022. URL https
11214
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2011666777779999 2019.
Jack Richter-Powee Yaron Lipman, and Ricky T. Q. Chen. Neural conservation laws: A divergence- free perspective. In S. Koyejo. S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 38075-38088, Curran Asso- ciates, Inc., 2022. https: /Rt.h/pp:e//ppreeednns....rrrree fi /ileeeeeeee
Luning Sun, Han Gao, Shaowu Pan, and Jian-Xur Xun Wang. Surrogate modeling for fluid flows based on constriied deep learning without simulation data. Computer Methods in Applied Mechanics and Engineering, 361:112732, 2020. ISSN 0045-7825, doi: https://doi.org/10..100
j.cma 12732.. URL https::///// www. www.nceiiecediieeoo rect. htw.o//nnnnerrrete.ooocceeeeeeee com/sscienne//rrtccle..
e11/8004578251930888
Fabricio Arend Torres, Marcello Massimo Negri, Marco Inversi, Jonathan Aellen, and Volker Roth. Lagrangian flow networks for conservation laws. In The Tweifth International Conference on Learn- ing Representations, 2024. URL /tt.httpee/////eeeerrreeee
wwwp//wwpwr.eeeeeeeeeeeeeeeoe.ew wet.i/ewwweeeeeeeeeeeeeee n forum www.roormmminnsssss-ssoss Nskkkkk
Alessandro Trenta, Andrea Cossu, and Davide Bacciu. Learning and transferring physical models through derivatives. Transactions on Machine Learning Research, 2026. ISSN 2835-8856. URL
https:://ppenreviee net forum? IIdrummmDdIII
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM J. Sci. Comput., (115):33355-3000 January 2021.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why PINNs fail to train: A neural tangent kernel perspective. J. Comput. Phys., 449 10768): January 2022.
Sifan Wang, Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. An expert's guide to training physics-informed neural networks, 2023. URL wwtpss://aaatttt. www./orggga...... abs/ 08468
Sifan Wang, Ananyae Kumar Bhartari, Bowen Li, and Paris Perdikaris. Gradient alignment in
physics-informed neural networks: A second-order optimization perspective, 2025. URL https:
7/arxiv. wwwp:///g//aaaaaaa.oo ablb///a222222 00604
Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non- adaptive and residual-based adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 403:115671, 2023. ISSN 0045-7825. doi: com/
Haixu Wu, Huakun Luo. Yuezhou Ma, Jianmin Wang, and Mingsheng Long. Ropinn: Re- gion optimized physics-informed neural networks. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (eds.), Advances in Neural Information Process- ing Systems, volume 37, pp. 10444-110532 Curran Associates, Inc., 2024. doi: 10.52202/ 079017-3508. URL ht / o
13

---

Published as a conference paper at ICLR 2026
Zixue Xiang, Wei Peng, Xu Liu, and Wen Yao. Self-adaptive loss balanced physics-informed neural networks. Neurocomput. 496 July 2022. ISSN 0925-2312. doi: 10.1016/j.neucom 2022.05.015. URL URLLtttpe///od...r//1000L.. 101.//1...00111111111 neucom. 2022. 05. 015
Yeonjong Shin Yeonjong Shin, Jérôme Darbon Jérôme Darbon, and George Em Karniadakis George Em Karniadakis. On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type PDEs. Commun. Comput. Phys., 28(5)20042-200 January 2020.
Xuan Zhang, Jacob Helwig, Yuchao Lin, Yaochen Xie, Cong Fu, Stephan Wojtowytsch, and Shuiwang Ji. Sinenet: Learning temporal dynamics in time-dependent partial differential equations. In The Twelfth International Conference on Learning Representations, 2024. URL https://
openreview net mett/Eorrum2dddSShEEhhi
Colby L Wight & Jia Zhao. Solving Allen-Cahn and Cahn-Hilliard equations using the adaptive physics informed neural networks. Commun. Comput. Phys., 29(3):93--954, June 2021.
Congcong Zhu, Xiaoyan Xu, Jiayue Han, and Jingrun Chen. Physics-informed temporal alignment for auto-regressive PDE foundation models. Forty-second Intermattionll Conference on Machine Learning, 2025. URL htt.htpp:: htt.::::///peerevieeee net net//ooruum2L==OOON
14

---

Published as a conference paper at ICLR 2026
A LLM USAGE
Large Language Models (LLMs) were used to polish and check the grammar of the paper, with very limited usage to rewrite some phrases for improved readability. No LLM was used to write entire parts of the text or for any other non-specified purposes. Every output of the LLMs was carefully checked to ensure its correctness and validity, adhering the LLM policy and Code of Ethics.
B MODEL DESCRIPTION
This Section provides further information CP models and their components, as well the employed baselines.
B.I PARTICULAR NCL APPLICATIONS
While PINNS have a straightforward implementation, NCL models require performing final trans- formation the network's s output. In particular, the MLP parametrizes antisymmetric matrix A, which depends on the input coordinates (4, 4r y). As described in Section 2 a row-wise divergence performed to obtain vector that has divergence with respect all the inputs. To obtain vector field that is divergence- free with respect to only a subset the input coordinates, it is sufficient consider smaller matrix A and apply the above procedure only for the interested coordinates. We provide some examples here.
In the case of the 2D time- dependennddnt Navier-Stoker equation, the incompressibility equation is given by
Ou. dx dy
Juy O,
V u (x.() (w)
where the field u depends on the three inputs (f, x, y). Here, the MLP parametrizes A as 2 x 2 antisymmetric matrix depending on (t, x. y). The row-wise divergence is taken only with respect to (x. y), so that the final output is divergence-fres in just those variables while depending freely on t. In the case the NS-Euler equations, the mass conservation for the density p is given by
Op Ōt
+ V (pr)) (i. (p. pm) 0.
(10)
To achieve this, the MLP parametrizes as 3 x antisymmetric matrix and the row-wise divergence taken with respect to all three variables (t,r,y). While divergence-free equation is guaranteed, other constraints can be added via losses similar to those of the PINN model. Finally, we consider the non-trivial case of the NCL model for the y- yevelocity equation in the Acoustics equations experiment in Section 3.3 which reads
Ou Ōt
Öp = 0. dy
(11)
To achieve this, it is sufficient for the MLP to output a 2 x 2 antisymmetric matrix (which is parametrized by only one output). Then, and p can be calculated as
Aoo
AAp
v = div (Ao.) = Øt
dy
(12)
8Aio
0Au
p div(r.)) (A1.)
0t
Öy
This way, the vector field (v, p) predicted by the NCL depends on all three variables (t, v), but is divergence-fin with respect only (t, v). While this example still simple write and involves variables, the above process can be generalized to any number and type of coordinates and for every possible subset of them.
B.2 SOBOLEV NORMS AND ALIGNMENT MODULES
In the mathematics literature, Sobolev spaces Maz ya 2011) the common choice to study and prove the existence and uniqueness of solutions to PDE problems (Evans) (022).. From a theoretical
point of view, it is not enough to measure distances in the 2 especially for convergence results. In
15

---

Published as a conference paper at ICLR 2026
addition, not all PDE problems admit smooth solutions, and notion of weak derivative is needed; hence, the necessity to consider a more restrictive space, that is, the Sobolev space W (S2) of functions u(x) with bounded norm
||((2))ll===mg = llx()))l +
(13)
where D'u(z) is the 1-th order differential of u. Weak derivatives of a function are defined such that they resemble integration by parts for every possible test function E Co (9t) with compact support in S2 that is
00 gos+aptan +. Js oz Ox2" D.... ..0
dx (-1)attmton [
go gorta2+akt..
dx (14) Ox an Ox o.... Ozh
Furthermore, while the derivatives and the boundary and initial conditions completely determine a function u, and are therefore sufficient to learn u (Trenta et al.] 2026), information on u is not sufficient to determine its derivatives, which fundamental for the dynamics and evolution of system. In practice, using the Sobolev norm to train a neural network has been shown to improve the performance, especially with low data availability (Czarnecki et al., 2017). In the context of physical systems, the derivatives are fundamental to ensure the physical system is learned and also to ensure the correct transfer of physical information between models (Trenta al.. (2026). For these reasons, the alignment module must contain the derivatives of the networks to ensure that the dynamics of the
system, predicted by the two networks, coherent. This explains why the alignment module with only u (OUTL) performs so poorly compared to DERL and SOB.
B.3 3 PDE RESIDUALS AND PINNS
In this Section, we provide more details on PDEs and their residuals. A PDE of order d is, in general, an expression that involves one or more partial derivatives, up to order function w:
Flu] F (u. Du, D2 D² ..., D²u) = O.
(15)
In the linear case, for example, it can be written as
Du
pPu Oua Oxr
Flu](2)) fxx + (o())u(()) + Σ ax(x)) Dz.
(x) +
(x)(
(x) +
(16)
PINNs implement the PDE itself as soft target in the loss. To do so, the PDE operator F, which takes as input a function and outputs another function, which we aim to be null, is calculated on the MLP itself. Derivatives of the MLP û, Dū,.. are calculated using Automatic Differentiation (Baydin et al. 2018), and the function F is calculated with them, obtaining the current PDE residual of the MLP itself
Flu) F (û, Dū, D'u....,D'i)
(17)
The norm of this residual is calculated and used to optimize the network. Ideally, if this loss term is 0 or close to it, the model itself should satisfy the PDE. This is also why we analyze it in Section (.44 see how well model is consistent with the underlying equation. Coupled with the loss for the initial and boundary conditions, this is the main component of Physics- Informed Neural Network (Raissi et al. 2019
B.4 GRADIENT- BASED REWEIGHTING
Gradient-bacee reweighting is a technique that automatically selects hyperparameters for the different loss terms based on the current gradient norm. We consider the approach described in Wang et al. (2023), where a model has parameters 0 and loss function
Lux (0) Aac Lic XBC + APDR (f).
(18)
16

---

Published as a conference paper at ICLR 2026
Then, the optimal weights calculated
drc Velccll +
VeCBCl
+ VPDpDE
V VoCicl
VaLic|| +
V VoC
+ V CpDE
ABc
(19)
VeCBcll
VoCiel +
V CBCll
+
V
ÂPDE
This formulation generalizes easily to multiple PDE terms and the alignment term.
B.5 ADAPTIVE POINT RESAMPLING (RAR)
For our adaptive point resampling, consider RAR (Wu al. 2023). After certain number of iterations, the PDE residuals of the system are calculated on grid of points. For each PDE, the k points with the highest residual are collected and added to the training data.
C IMPLEMENTATION DETAILS
In this section, we provide details the practical implementation of the CP model, from the discretization of the losses the computational times and the tuning process.
C.1 LOsS DISCRETIZATION
Since not possible to calculate the true T.2 losses and distances for the models, we discretize them as is usually done in the literature. In general, assume have 4 datasets:
The first contains ND points in the interior of the spatio-temporal domain [0, T] x Ω with no given supervised targets:
Do = {(t,, 2.), E [0, E Ωi = Np).
(20)
Points in this dataset are often called collocation points in literature and are used in the PDE and alignment loss. The intuition is that, if these losses are close to zero, the model satisfies the PDE (or the networks are aligned) for sufficient points in the interior of the domain. The PDE residual loss is then approximated as
| Mi [F
[F(ū)(4, a.)1?
(21)
ND
while the alignment loss components discretized in a similar way
No c (1:, xx)) - a2:(1;
h u2 22212 221122(0.77x02))
R
ND
(22)
1
D - D'
j (t, x.) D xt))11
ND
The second contains N, points in the domain Ω at time t and is used to enforce the initial condition u(0, x) = g(x):
Dj = {( (0, x.). æ; E g2 e R", 1.... Nr
(23)
The IC loss is then discretized as
N
a - 11
2
x;;) - g(xx)))1²
(24)
Nr
17

---

Published as a conference paper at ICLR 2026
Table 4: Time in seconds required for one step of training.
Model
Taylor-Green vortex
Kovasznay flow
Acoustics
Euler gas
MaghetotHydroDyymmiss
PINN
0.0420
0.0480
0.0210 0.0618
0.2097
0.1093
NCL
0.1807
0.1174
0.6889
0.2553
2xPINN (OUTL)
0.1034
0.0840
0.2884
2xPINN (SOB)
0.1069
0.0947
0.3281
0.1657
2xPINN (DERL.)
0.1039
0.0933
0.3101
0.1611
PI+NCL (OUTL)
0.1431
0.1562
0.5795
PI+NCL (SOB)
0.1483
0.1543
0.6584
PI+NCL (DERL)
0.1413
0.1626
0.6563
2xNCL (SOB)
0.6896
2xNCL (DERL)
0.6858
3xPINN (OUTL)
0.0855
0,4366
3xPINN (SOB)
0.0898
0.4366
0.2025
3xPINN (DERL)
0.0888
0.4120
0.2002
3xNCL (OUTL)
0.1222
3xNCL (SOR)
0.0621
3xNCL (DERL.
0.0699
4xPINN (SOB)
0.2341
4xPINN (DERL)
0.2207
The third one contains NE points on the boundary 02 and for different times [0, [0, It used to enforce the boundary conditions B/u]] x) = b(t,
Dj = ((t,, z;), b(t,, a))), t; E [0, E Ω b(t,()))))) E
NB). (25)
Boundary conditions can be of many forms, such as Dirichlet BC, which impose u(t. x)) = b(t, æ), or Neumann BC, which impose condition on the derivative (t, b(t, where ñ is the normal vector at the boundary. The loss, in this example for Dirichlet conditions, is then discretized as
N NB Σ
B[a] - (1(2a(()))
x;) b(t,, xx))
(26)
The last one, used only evaluation, contains a set of points the spacetime domain (1;, ax)) and the corresponding true solution of the system u(t;, m:):
Ds {{ (f:, R,),, u(t;, 2.)),,t & E [0, 7). x. E Ω u(t;a;) R", 1..., Ns). (27)
This dataset is used to compare the final solution of the models against a true target, never seen during training. In this case, the loss is discretized as usual with
Ns ||u(ts, x.) - a(t;, 2;)1
||u - 0012A(0,27x08)
Σ
(28)
Ns
C.2 COMPUTATIONAL TIME AND RESOURCES
In Table we provide comparison of the computational times of the different methods on all experiments. We reported the time to perform one epoch of training for each experiment. All models were trained on a single NVIDIA H100 GPU with 80GB of memory.
C.3 MODEL ARCHITECTURES
We report the architectures of the models employed in Table
D ADDITIONAL RESULTS
In this Appendix, we provide further experimental details and results for the experiments Sections
3 and 4
18

---

Published as a conference paper at ICLR 2026
Table
MLP architectures used in the experiment.
Model
Taylor-Green vortex
Kovasznay flow
Acoustics
Euler gas
MighetooHydnoDyaamiiss
PINN
Layers and units
2x64
4x64 Tanh
4x64 Tanh
8x512 Tanh
8x256
Activation
Tanh
Tanh
NCL.
Layers and units Activation
2x64
4x64 SoftMax
4x64 SoftMax
8x512
8x256
SoftMax
SoftMax
SoftMax
CP-PINN modules
Layers and units Activation
2x64
4x64 Tanh
4x64 Tanh
8x256 Tanh
8x256
Tanh
Tanh
CP-NCL modules Layers and units Activation
2x64
4x64 SoftMax
4x64 SoftMax
8x256 SoftMax
8x256
SoftMax
SoftMax
ONZITAN Gredd
GN
35
34
4000
34
COR
(b) t = 2.
FRN
FRIU--A
ICL
(c) = 4.
389
13
(d)
Figure 6: Taylor-Green vortex. Model pointwise enror for the solution (u,p). Only the best alignment for each module choice shown.
D.1 TAYLOR- -RRENN VORTEX
The analytical solution of the Taylor-Green vortex (Chorin 1968) used in experiment is given by
(1x () = [sin x cos y, cos x sin y) e 2et p(1, y, 1) 1 (cos(2xr) + cos(2y))e 4rt (29)
Figure 6 shows the model errors at different times.
19

---

Published as a conference paper at ICLR 2026
140
340
327
Figure 7: Kovasznay flow experiment. Prediction errors for PINN (and its variations), NCL, and CP models (one for each combination of modules). The lowest errors are in blue.
Laver
101
-0.50000 025 0002 0.005
0.015
0900000-00
0.000 0.225
0550 C.055
0.550-0.025 0.09
6,023 0,000 0.115
18.00-00 025
3000
0.555
0.050
0375
mEn
(a) Kovasznay flow 2xPINN (DERL)
Layer
Loyer
Layer 2
Layer
-0.24
11:00
2112
11.04
-8.08
-0.22
1105
11.02
2118
0.00
0122
6164
-0UR
-0.0
0.50
0.17
013
(b) Kovasznay flow PINN
Figure 8: Gradient histograms. Each plot contains the histograms for the distribution of the gradients propagated at each layer of CP and PINN at the beginning of training, similarly to Wang et al. (22TT)
D.2 KOVASZNAY FLOW
The analytical solution to the Kovasznay flow (Drazin & Riley 2009) is given by
X Ae sin sing(
u(r,y) 1 e Ae cos cos())))) where
p(a, v)
2200)
(30)
2n
To + 2 this case, the vorticity given by Se he sin sin()
Figure
20
2x
shows the errors of the employed models in the domain.
Gradient Analysis. Here, we provide a similar analysis to the one Section 3.4for for the Kovasznay flow experiment, comprising 2 equations and modules. Figures Sal and 8b show the distribution of the propagated gradients in a training step at the different layers of the MLP, respectively for the CP-2xPINT (each row is module) and PINN models. The gradients for each loss term represented with a different color. We clearly see the effect learning the system in a compositional manner: the distribution histograms of the gradients of the loss terms are much more aligned than in the PINN model. This means that each loss can be learned more effectively without additional regularization (Wang al. 202T
D.3 ACOUSTICS EQUATIONS
The initial condition of the Acoustics experiment in Section 3.3 is given by
p(0 x. y) = 1 1 cos
0.2 (v y2
IV
x2
+ 2
0.5| -
< 2,
(31)
20

---

Published as a conference paper at ICLR 2026
PINN
PINN-Grad
NCL
3XPINN (DERL+Grad)
CP (DERL (DERL+GTad)
0.30
055
5
0.0
64
-0%
0.0
01
-0%
50
05
-0
05
05
-51
00
04
(a) 0.
PIN
PINNeGrad
NCL
CP-3XPNN (DERL+Gram) GGrad)
CP-3xNCC (DERL+GraD
0.0
0.0
0.0
0d5
(b) 0.08.
(DERL+Grad)
055
0.20
0.0
0.5
-0.3
0.0
63
-0.5
00
0.
-0.
00
0.5
5
0.
05
-0.3
0.0
93
(c) = 0.16.
PINN
PNN--Gra
NCL
CP CP33XPNN IDeRL+
CPOCNNC (DERL+Srad)
02
-05
-0.
0.0
-0.5
-5.5
-0.5
0.
-0.5
0.0
0.5
0.24.
Figure 9: Acoustics equations. Model pointwise error for the solution (p, f). Only the best alignment for each module choice is shown.
and null velocities in the whole square. Figure SourreSshe shows the model errors at different times.
D.4 NS-EULER GAS EQUATIONS
provide further details on the setup of the experiment in Section The initial conditions are given by
Po(x, y) (sin(2mx:) + sinn)) 1,
uo(x,
[ c cin((2m))
/2
(32)
To embed points from the 2D unit square [0, 1]2 to the two-dimensional torus T, the following function is used:
:: [0, 1]2 - T i(z, y) (cos(2 sin(2:..)), cos(2rr)) (2n(2))))
(33)
which also automatically imposes periodic boundary conditions on the square.
Figure 10 shows the model emror different times. We that the PINN and NCL models suffer in certain regions of the domain, while our models do not. this case, Gradient- based reweighting worsens the performance of PINNs.
21

---

Published as a conference paper at ICLR 2026
86
2.8
2.8
84
244
27
25
1,0 00
00
05
1,0
00
0.5
1.0 00
05
3.0 60
1.0 00
03
3.0
(a) = 0.
PNNN
NCL
PKN--NCL NCL (DERLI
2xPANN (DERL)
2XNOL (DEALI
3UPANN (OERL)
20
3.6
35
80
(b) 0.1.
5.4
PINN+-rad
PINI-NCL (DERLI
2XRNIN EXNCL
2.8
2.0
2.4
27
05
1005
32
00
05
10
35
L602
0.
10.00
65
20:00
33
100
(c) t = 0.2.
PMN
NCL
PINK-NCL (DERU)
2*PINN IDEAL!
EXNCL (DERL)
2xPNN (DERL)
24
38
80
(d) 0.3.
Figure 10: NS-Euler gas equations. Model pointwise error for the solution (p, u. p). Only the best alignment for each module choice is shown.
22

---

Published as a conference paper at ICLR 2026
(a) t= 0.1.
Eex
(b) 0.3.
(c) L = 0.5.
Figure 11: MHD equations. Model pointwise error for the solution.
D.5 MACNETOHYDRODDNAANTN
Here, we provide additional details on the experiment in Section 4.2
The initial condition is given by
7
p(0, x, 4m
u00,z, x. |- sin (2Ty), sin (2Tx))
P(O, x. v)
(34)
4T
cos( 4TT
cos (2s))
A(0, x, 4TH 4m
2T 4m
B(0, x, y) = V x A
where A(t, x. is the magnetic potential Gruber & Rappaz 1985). Similar to the NS-Euler experiment in Section 4.1 to satisfy the periodic boundary condition, we embed the square domain [0, [] x [0, in the two-dimensional torus T with equation 33 from Appendix D.4 The data generation done as in Giopakumar et al. (2025)
Figure 11 shows the prediction enrors of the models in the domain different times in [0, 0.5].
FURTHER ANALYSIS ON THE ALIGNMENT MECHANISM
Loss curves. We start by showing how the alignment loss, the prediction error, and the PDE residuals for transferred constraints are related to each other in practical examples. Figure 12 reports smoothed loss curves for the CP models on the Acoustic equation experiment of Section 3.3 These losses calculated the inference module, which learns equation (A.P) in table but never learns equations (A. Vx) and (A. Vy) directly. These PDEs are learned by the other two modules. The inference module leams satisfy these constraints only by aligning with the other two modules. Figure 1 shows smoothed loss curves for the CP models on the Kovasznay flow experiment in similar manner.
We now investigate how well the modules are aligned during training by keeping track of their prediction errors for the common variables during the whole training process. This shows how
23

---

Published as a conference paper at ICLR 2026
MP-3XPINN (DERL) Losses
MP-3XNCL (DERL) Losses
10-
(A. residual
(A. residual
(A. residual
(A. residual
Alignment Loss
Alignment Loss
Prediction Error
Prediction on Error
1
Loss
10-3
20000 40000 60000 80000 100000 Step
20000 40000 60000 80000 100000
Step
Figure 12: Smoothed loss curves of CP models on the Acoustic equation experiments. Losses are calculated on the inference module, which never learns equations (A.Vx, and (A. Vy) directly.
MP-2XPINN (DEKL) ) Losses
MP-PINN NCL (DERL Losses
(KEM) residual
(KEM) residual
10a
(KFI) residual
(KF-) residucl
Alignment
Alignment
10
Prediction Eiror
Prediction Eiror
10
10
10
20000
4UX0X
600000
BU0O0
1000000
200000
400000
600000
80000
Stap
Step
Figure 13: Smoothed loss curves of CP models on the Kovasznay flow experiment. Losses are calculated on the inference module, which never leams equation (KFII) directly.
24

---

Published as a conference paper at ICLR 2026
Prediction losses for MP-2xPINN GEaL)
Prediction losses for MP-PINN NCL (DERL+ OBERL+GraDe Gradi
Momenlum
10'
Momenium mocule
Divergance module
Divergence module
10-3
10ah
1.0
10-
10
10.5
10-4
10a
10
200c0
000
50000
80C.0
190000
20000
4C0C0
60000
80000
100000
Tr Training
Training
(a) CP-2xPINN
(b) CP-PINN+NCL
Figure 14: Prediction errors on the shared variables for the CP modules on the Kovasznay flow experiment
well their predictions are aligned during the whole process. If the alignment mechanism is working as desired, the loss curves should be very similar, meaning that the modules working learning together. For the Kovasznay flow experiment, we report the errors for u for the Momentum and divergence modules in Figure [4] for both the CP-2xPINN and CP-PINN+NCL models with DERL alignment and gradient-basee reweighting. For the Acoustic equations experiment, we have that the pressure and velocity x modules share variables (p. u), while the pressure and velocity y modules share (p, v). Hence, we have three plots: one shows the prediction errors for all three modules on p, one shows the prediction errors for the first and second module on (p, u), and the third shows the precision errors for the first and third module on (p. v). These results are available in Figure 15 for both the CP-3xPINN and CP-3xNCL, models with DERL alignment and gradient- based reweighting.
and
These plots clearly indicate that the losses are very similar throughout the whole training process, indicating that the alignment mechanism is making the modules collaborate closely to obtain the final solution.
Measuring discrepancies. We now show from a numerical point of view how close the functions predicted by different modules are in practical terms. We perform an additional test on the Kovasznay Flow and Taylor-Green experiments from Section 3 At the end of training, we predict the solution for the variable which is a vector containing the and y components of the fluid velocity, shared among the two modules of the CP models. We measure the following quantities:
The distance between the two modules, normalized with respect to the true solution norm MMnatate 1 - 222 this indicates how different are the two predicted solutions, scaled to the M magnitude of the true function.
For a point-wise comparison, we consider the symmetrized Percentage Error: for each point 1(1,a) -Vmaai:
in the domain, we calculate [[ (ff.z]]| [[f.ac] which measures how far are the two predictions compared to each other. We then calculate the mean and maximum value.
We perform this test for both the CP-2XPINN and CP-PINN+NCL models. Results are available in Table From the scale of the metrics, we conclude that the difference between the two modules is negligible compared to the scale of the solution with respect to each other. We have not observed any cases where the modules disagree with each other, as they are practically the same function when it comes predicting the solution.
25

---

Published as a conference paper at ICLR 2026
Convergence of Losses for MP-3XPINN (DERL+Grad)
Error on (p, u)
Error on (o. v)
Error on p
Pressure moedule Velocity x madule
Pressure module Velocity y module
Pressure module
Velocity x module
10 1
Velocity y module
10-2
10-3
50000 Training step
100000
0
50000
100000
0
50000 Trining step
100000
Training step
(a) CP-3xPINN
Convergence Losses for MP-3XNCL (DERL Grad)
Error on (p. u)
Error on (p. v)
Error on
Pressure moedule
Pressure module
Pressure module
Velocity x madule
Velacity y module
Velacity x medule
10-1
Velocity y module
Loss
10
10-3
50000
100000
50000
100000
50000
100000
Training step
Training step
Training step
(b) CP-3xNCL
Figure 15: Prediction errors the shared variables for the CP modules on the Acoustic equation experiment
Table 6: Normalized 2 distance and symmetrized Percentage Error metrics for the modules on predicting the Kovasznay flow and Taylor Green vortex experiments.
Kovasznay Flow Experiment
Model
Normalized L distance
Mean sPE
Max sPE
x 10
x10-7
x10
CP-2XPINN (DERL+Grad)
2.758
6.095
2.954
CP-PINN+NCL (DERL+Grad)
2.210
4.874
1.729
Taylor-Green vortex Experiment
Model
Normalized T.2 distance
Mean sPE
Max sPE
x10 4
x10-7
x10-6
CP-2XPINN (DERL+Grad)
8.910
2.762
2.102
Co-PIN-NNCLDEEL--Grrdd)
10.47
3.187
4.693
26

---

Published as a conference paper at ICLR 2026
Table 7: Results for different values of Aaligo on two CP models.
Acoustic equations
CP-3xNCL (DERL)
L error
max error
x10-3
x x10-1
aaiggn 1
5,165
1.547
Aalign 0.1
3.164
1.257
Aalign 0.01
2.460
1.089
Aalign 10.
12.46
4.388
Table 8: Results for the Acoustic equation experiment where NCL and PINN have the same parameter budget of CP during training. In these cases, PINN and NCL have times the number of parameters of CP during inference.
Acoustic equations
Model
I error
max_err
x10-3
x10-1
PINN
5.671
1,462
NCL
2.940
1.169
CP-3xNCL (DERL+Grad)
2.718
1.121
F ABLATIONS
Alignment coefficient Aatign- We experimented with different values of Aalign to analyze the sensi- tivity of CP with respect to this hyperparameter. We test values of Aalign between 0.01 and 10. Results for the Acoustics experiment (CP-3xNCL model with DERL alignment) are available in table
As we can see, Aalggn regulates the importance of the alignment between models. While its influence depends on the specific task and modules, the Aalign 1 works well on average in all benchmarks. Highly specialised tuning of the hyperparameter can lead to modest improvement, as in the acoustic equation case above.
Number of parameters in the baselines. Since the CP model employs more than one MLP during training, the total number of parameters during this phase is increased. At inference, where the models can be used indefinitely, CP uses a single module, and hence will have the same number of parameters and expressiveness as an NCL PINN module. In this ablation, we increase the number of units in the NCL and PINN layers to match the parameter budget of the CP model during training. We experimented with the Acoustics equations and the CP model with 3 NCL modules. We remark that, during inference, this means that PINN and NCL. will use 3 times the parameters of CP. Results are available in Table g As we can see, our models perform better even though they use fewer parameters during inference, further showing the importance of modularization in these applications.
Statistical Significance. To show robustness with respect to different initialization and stochasticity in optimization, we perform additional runs with different seeds on the Acoustic equations experiment. We experiment with the PINN and NCL baselines as well as the CP-3xNCL model, DERL alignment, and gradient-based reweighting. Results averaged over three seeds are available in Table
As we can see, the results clearly indicate that our CP model with NCL modules still the best model by far, with lower L and maximum errors in the domain. Furthermore, the standard deviations are smaller, indicating more robustness to different initializations.
G FURTHER THEORETICAL DISCUSSION
In this Section, we provide additional theoretical understandings and intuitive insights on why our modular approach is easier optimize compared to PINN model.
27

---

Published as a conference paper at ICLR 2026
Table 9: Results for the Acoustic equation experiments with different seeds. Results are averaged over 3 runs. We report mean and standard deviations.
Acoustic equations
Model
L error
max err
x x10
x10
PINN
5.090
2.247:10.418
NCL
3.540=0.431
1.610=0.082
CP-3xNCL (DERL)
1..73
On the difference between PDE residual terms and alignment losses PDE residual term imposes constraints that link together the partial derivatives of MLP. For example, PDE such Oy A Vu 0 imposes that the time derivative of the network must be equal to a linear combination of the spatial ones. As observed in Krishnapriyan et al. (202T), this type of constraint is often ill-conditioned and can lead to numerical instabilities. The same work also shows that these PDE terms produce very complex and non-smooth loss landscape. Hence, our idea is avoid combining multiple PDE losses in a single MLP. ComPhy tackles subproblems (simpler than the whole system) in each model, and then transfers information across models through alignment.
as ăr
The alignment loss, on the other hand, does not link together the partial derivatives or the outputs of MLP. Instead, uses the derivatives or outputs of one module as a target for those of another module in supervised manner (to be precise, this is similar to a distillation between the two models). Hence, we expect this type of loss to reduce the complexity of the problem and to simplify the optimization process.
Physical constraints and properly aligned modules. When two models are properly aligned, that the functions predicted by the modules
is, when the SOB DERL alignment loss is close are very close in the Sobolev W112 space (Trenta et al. 2026). As an example, we use the one in Section 2.3 on the Navier-Stokes equations. After training. module 1 has learned the momentum equation, while module 2 has learned the divergence equation. Since Umodule & Umodule and D 1 Dumddall 2, we calculate the divergence of module 1, we can say that -mmoulle & V Umodule ~ 0. This means that correct alignment allows module to satisfy the divergence equation because module does so. This intuitively provides reason why module is still capable of learning the system while seeing only one equation. We experimentally validated this claim in Appendix B
28


Reviews:
[{"review_id": "L8fe3MWOgP", "reviewer": "Reviewer_amcN", "summary": "This paper introduces ComPhy (CP), a modular framework for solving systems of Partial Differential Equations (PDEs) using machine learning. The core innovation involves assigning each PDE in a system to a dedicated learning module (such as PINNs or Neural Conservation Laws) and connecting these modules through an end-to-end alignment mechanism. The alignment process enforces consistency between modules that predict the same physical variables, with particular emphasis on derivative-based alignment losses (Sobolev norm and derivative-only alignment). The authors demonstrate that this compositional approach outperforms standard methods where a single model learns all PDEs simultaneously. Experiments span systems ranging from two to five equations, including Navier-Stokes, acoustic wave equations, and magnetohydrodynamics, showing consistent improvements in accuracy. The paper also provides gradient analysis suggesting that CP's modular structure leads to more balanced gradient distributions during training, which may explain its superior performance.", "questions": "Statistical significance and error bars\nThe results tables (Tables 2 and 3) report point estimates without error bars or confidence intervals. Given that neural network training involves stochastic elements (random initialization, batch sampling), reporting means and standard deviations across multiple runs would strengthen the claims.  \n\nAblation on alignment losses\nThe authors should include ablation studies showing performance across a range of λ_align values, demonstrate the method's robustness (or lack thereof) to hyperparameter choices\n\nComputational Efficiency and Cost–Performance Analysis \nThe training time overhead of CP models is mentioned but not quantitatively justified. It would be useful to discuss whether the additional training cost is proportionate to the observed performance improvement. Presenting these metrics side by side would improve the completeness and transparency of the experimental analysis.\n\nParameter-matched baseline  \nCompare against a single PINN baseline that is parameter-matched to the full CP system (same total parameter count). Alternatively, show scaling curves where single PINN and CP models are compared across increasing total parameter budgets. This will clarify whether gains are due to modularity or simply model capacity.\n\nAlignment learning\nRecent studies [1][2] have incorporated alignment learning into PINN-like frameworks, demonstrating its potential to enhance physical consistency and optimization efficiency. Therefore, the authors should discuss how their proposed method relates to these works, highlighting the key differences or advantages.\n\n[1]Gradient Alignment in Physics-informed Neural Networks: A Second-Order Optimization Perspective\n\n[2]Physics-informed Temporal Alignment for Auto-regressive PDE Foundation Models", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Strengths\n\nThe paper demonstrates several notable strengths that support its contribution to the community.\nClear motivation and intuitive approach: The paper effectively motivates the problem of solving coupled PDE systems and presents an intuitive solution. The compositional structure mirrors the mathematical structure of the underlying physical system, making the approach both theoretically appealing and practically sensible. The gradual build-up from methodology to the concrete Navier-Stokes example in Section 2.3 aids understanding.\nStrong and consistent empirical results: The experimental evaluation is comprehensive, covering multiple physical systems with increasing complexity. \n\nThorough experimental methodology: The authors compare against multiple relevant baselines including PINN with gradient reweighting and adaptive point resampling. The experiments are well-documented with detailed problem setups, boundary conditions, and reference solutions in the appendices. \n\nValuable gradient analysis: Section 3.4's gradient histogram analysis provides meaningful insight into why the modular approach succeeds. The observation that CP produces more balanced gradient distributions across layers compared to standard PINNs offers an empirical explanation for the performance gains and could inform future research.\nGeneralization beyond divergence-free equations: Unlike NCL which is specifically designed for divergence-free fields, ComPhy's framework applies to general PDE systems, demonstrating particular value in experiments like acoustics and MHD where NCL-only approaches would be insufficient.", "weaknesses": "Weaknesses and Concerns\n\n Insufficient analysis of hyperparameter selection and sensitivity\nThe paper does not provide clear guidance on choosing the critical hyperparameter λ_align. While Table 2 and Table 3 show results with fixed hyperparameters, there is no ablation study examining sensitivity to this choice or methodology for setting it. \n\n Module assignment strategy not systematically addressed\nThe paper does not provide principled guidance on how to assign PDEs to modules. For the NS-Euler experiment (Section 4.1), the authors test multiple configurations (2xPINN, PINN+NCL, 2xNCL, 3xPINN) but offer no systematic approach for making this choice. Different assignments can lead to different architectures, but the selection process appears ad-hoc.\nThe paper would benefit from either developing heuristics for module assignment (e.g., based on PDE type, coupling strength, or variable sharing patterns) or demonstrating that the method is robust across reasonable assignment choices.\n\n Comparison fairness and architecture choices\nThe baseline single PINN appears to use the same architecture size as individual CP modules (Table 5), meaning the total CP model has substantially more parameters across all modules. It is unclear whether a larger single PINN with comparable total parameter count would close the performance gap. I wonder whether the observed gains arise specifically from the modular structure or merely from the increased model capacity.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "XcdLvfeOCa", "reviewer": "Reviewer_AyTi", "summary": "This paper presents ComPhy, a novel modular framework designed to leverage the inherent physical structure of the problem to solve systems of PDEs. ComPhy assigns each equation a dedicated learning module, like PINN or NCL, and introduces an end-to-end alignment mechanism to enable knowledge transfer between different modules. The results show that it outperforms state-of-the-art approaches where a single model is trained on all PDEs at once.", "questions": "See weaknesses.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- From my perspective, this is a novel method designed for solving systems of PDEs. It is an interesting research field which have not been explored.\n- The proposed method is novel and seems elegant for solving systems of PDEs. The experiment results also demonstrate that it outperforms plain PINNs.\n- The paper is well written and easy to follow.", "weaknesses": "- It may be hard, but some theoretical understandings, even intuitive ones, could make the method more convincing.\n- After training, we can use a subset of the trained networks to predict all physical variables. Can some networks sharing the same variables have conflicts?\n- The paper lacks analysis on the efficiency of the model. For a system of N PDEs, each network requires (N+2) loss terms, and the whole system requires N*(N+2) loss terms; how does it influence the training time compared to plain PINNs.\n- Can we replace the current PINNs and NCLs with neural operators? The current framework does not seem to support such modules.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "v7bL9HbhIy", "reviewer": "Reviewer_Jj5v", "summary": "The paper introduces ComPhy, a framework for solving systems of PDEs. Instead of training a single monolithic model on all equations, ComPhy assigns each PDE to a dedicated learning module (e.g., a PINN or NCL) and enforces alignment losses between modules that share physical variables. These alignment losses encourage consistency across modules and improve convergence.", "questions": "1. Managing multiple interacting modules may increase computational and memory overhead, particularly for systems with many PDEs. It would be helpful if the authors could clarify how they address this issue.\n\n2. Since the overall objective combines both alignment losses and module-specific PDE/BC/IC losses, it would be useful to report how sensitive the method is to the relative weighting of these terms. Does performance degrade significantly if the alignment coefficient is varied?\n\n3. Can the modular design generalize to systems where PDEs share only partial or implicit variables?\n\n4. The figures, particularly Figure 1, could be made more intuitive and easier to interpret.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The paper introduces a novel modular design. The decomposition of multi-PDE systems into specialized modules is elegant and well-motivated both computationally and physically.\n\n2. Introducing Sobolev-inspired alignment effectively transfers physical information between modules and leads to empirical gains.\n\n3. Gradient distribution studies convincingly explain why ComPhy’s modular approach stabilizes training compared to conventional PINNs.", "weaknesses": "See questions below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "RS05mBq7VG", "reviewer": "Reviewer_M37D", "summary": "This paper presents a modular framework to tackle the optimization problem of PINNs caused by multiple loss constraints, which is named ComPhy. Specifically, ComPhy proposes to split multiple loss functions into several subsets with shared physical quantities. Then, several modules are configured as PINN or NCL models, and each module will be optimized based on one subset of loss functions. Besides, an alignment loss is newly proposed to align the shared physical quantities in different modules. As for inference, ComPhy only needs to infer the module with complete quantities once. Such a modular framework is expected to ease the difficulty in joint optimization of multiple losses and ensure the physical alignment in the final results. The authors provide sufficient experiments to verify the effectiveness of the proposed ComPhy.", "questions": "Please see Weaknesses. To highlight, I think the authors should answer the following questions carefully:\n\n-\tWhy does ComPhy work well?\n\n-\tHow to decide the configuration of ComPhy?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 1, "presentation": 1, "contribution": 2, "strengths": "(1)\tI think the proposed idea is interesting and novel, especially in assigning different loss subsets to different modules.\n\n(2)\tThe experiments are sufficient to deliver a comprehensive evaluation of the proposed method.\n\n(3)\tRich implementation details are included.", "weaknesses": "### (1) The motivation of ComPhy. Why does this design work?\n\nI think the authors fail to elaborate on why ComPhy performs well. The only statement is in Lines 99-101, that is, “State-of-the-art models like PINNs may suffer from optimization problems when multiple PDEs are involved. ComPhy avoids this issue by using different modules to optimize the different PDEs of the system separately.” However, suppose one of the modules in ComPhy contains complete physics quantities. In that case, it will also be optimized by the newly added alignment losses, which is still a multiple-loss optimization problem. Why is the module optimized in this way better than the model optimized from multiple PDE losses? I think an intuitive understanding is required.\n\nBesides, let us consider the example in Section 2.3. If ComPhy employs two PINN models as modules, at the beginning of training, it is really hard for the second module optimized with Eq.~(8) to generate a reliable solution. Why can the alignment loss help the first module be optimized better?\n\nTherefore, I cannot understand why ComPhy can help with the training. **Maybe the visualization of training curves (training, alignment and test losses) can be a good choice for elaboration. In addition, some theoretical analyses or thought experiments are expected.**\n\n### (2) Too many unjustified or vague claims.\n\nI think this paper contains many unsupported claims or statements, which seriously damage the scientific rigor. Here are some examples:\n\n-\tAbstract: “CP is the first approach specifically designed to tackle systems of PDEs”. Suppose the authors refer to “systems of PDEs” as the combination of multiple PDE equations. In that case, I think there are many related works that tackle the optimization problem of balancing multiple PINN losses, such as [1].\n\n-\tLine 99: “State-of-the-art models like PINNs may suffer from optimization problems when multiple PDEs are involved”. What kind of “optimization problems” do you mean? I think if the authors cannot detail this statement, the motivation of this paper is unclear.\n\n-\tLine 111: “The solution to a PDE system is unique only if all the PDEs are satisfied at once (Evans, 2022).” Although this is not a core statement, it should be noted that many PDEs contain multiple solutions.\n\n-\tEq. (3): The authors do not provide a clear definition for the L_{align}, since in Eq. (3), each row defines one type of L_{align}. I do not know what the final version used in Eq. (4) is.\n\n-\tTable 1: It is really hard to understand the last column. After a long time thinking, I understand that each row of the last column represents one configuration in ComPhy. I think a more direct description is required. \n\n-\tAll the numbers, like 100.000 or 600.000, should be 100,000 or 600,000.\n\n-\tLine 356: “The authors show that when PINNs are optimized correctly, the gradients tend to be evenly distributed across all layers.” This claim is not correct, since in Wang et al. (2021), the main focus is on the imbalance among different losses. Thus, this statement should be “the gradients of multiple PDEs”.\n\n[1] Wang et al. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 2022.\n\n### (3) How to decide the configuration of ComPhy, such as pure PINNs, pure NCLs, or a combination of PINN and NCL? \n\nAs presented in Table 2, different configurations lead to quite different results. When using ComPhy, it may take a long time to tune the concrete configuration of the modular framework.\n\n###  (4) About related work. \n\nI think this paper is not related to the neural operator, which is purely data-driven. The authors should spend more time reviewing papers about PINN optimization, such as [1,2,3].\n\n[1] Wang et al. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 2022.\n\n[2] Daw et al. Mitigating propagation failures in physics-informed neural networks using retain-resample-release (r3) sampling, ICML 2023.\n\n[3] Wu et al. RoPINN: Region Optimized Physics-Informed Neural Networks, NeurIPS 2024.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)