GenCP: Towards Generative Modeling Paradigm of Coupled physics
Authors: Tianrun Gao
Haoren Zheng
Wenhao Deng
Haodong Feng
Tao Zhang
Ruiqi Feng
Qianyi Chen
Tailin Wu
Venue: ICLR 2026 Poster (2026)
Keywords: Coupled Physics Simulation
Flow Matching
Operator Splitting
Joint Sampling
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: GenCP introduces a generative paradigm for coupled physics simulation that not only offers excellent theoretical guarantee, but also provides concrete guidance for practical implementation.

Abstract: Real-world physical systems are inherently complex, often involving the coupling of multiple physics, making their simulation both highly valuable and challenging. Many mainstream approaches face challenges when dealing with decoupled data. Besides, they also suffer from low efficiency and fidelity in strongly coupled spatio-temporal physical systems. Here we propose GenCP, a novel and elegant generative paradigm for coupled multiphysics simulation. By formulating coupled-physics modeling as a probability modeling problem, our key innovation is to integrate probability density evolution in generative modeling with iterative multiphysics coupling, thereby enabling training on data from decoupled simulation and inferring coupled physics during sampling. We also utilize operator-splitting theory in the space of probability evolution to establish error controllability guarantees for this “conditional-to-joint” sampling scheme. We evaluate our paradigm on a synthetic setting and three challenging multiphysics scenarios to demonstrate both principled insight and superior application performance of GenCP. Code is available at this repo: https://github.com/AI4Science-WestlakeU/GenCP.

Published as a conference paper at ICLR 2026
GenCP: TOWARDS GENERATIVE
MODELING
PARADIGM OF COUPLED PHYSICS
Tianrun G2 Haoren Zheng" Wenhao Deng"! Haodong Feng'
Tao Zhang 1,4 Ruiqi Feng' Qianyi Chen Tailin Wu
'Department of Artificial Intelligence, School Engineering, Westlake University; *Department of Geotechnical Engineering, Tongji University;
Global College, Shanghai Jiao Tong University;
*State Key Laboratory of Advanced Nuclear Energy Technology, Nuclear Power Institute of China; [gactianrun, chenqianyi, wutailinjewest.eee edd.. cn
ABSTRACT
Real-world physical systems inherently complex, often involving the coupling of multiple physics, making their simulation both highly valuable and challeng- ing. Many mainstream approaches face challenges when dealing with decoupled data. Besides, they also suffer from low efficiency and fidelity in strongly cou- pled spatio-temporal physical systems. Here we propose GenCP, a novel and elegant generative paradigm for coupled multiphysics simulation. By formulat- ing coupled-physics modeling as probability modeling problem, our key inno- vation is to integrate probability density evolution in generative modeling with iterative multiphysics coupling, thereby enabling training on data from decou- pled simulation and inferring coupled physics during sampling. We also uti- lize operator-splitingg theory in the space of probability evolution to establish sampling scheme.
error controllability guarantees for this "conditionallto-jonn We evaluate our paradigm on a synthetic setting and three challenging multi- physics scenarios demonstrate both principled insight and superior application performance of GenCP. Code is available at this repo: github.oom/Al4Scienne WestlakeU/GenCP
INTRODUCTION
Most real-world physical systems are governed by the intricate interplay of multiple physical pro- cesses spanning diverse disciplines (Yang et al., 2025). As a result, multiphysics problems are widely recognized as both fundamental and practically valuable. They arise across a broad range of applications, including aerospace engineering (Malikov et al. 2024 Hu al., 2024). biological engineering (Pramanik et al. [024] Gerdroodbary & Salavatidezfouii 2025), and civil engineering (Wang et al., 2025] Wijesooriya al. 22220). Despite their importance, accurate simulation cou- pled physics remains notoriously difficult (McCabe et al.) 2024)., largely due to strong cross-physics interactions and the resulting high system complexity, which are far more challenging to model than in single-physics settings.
Numerical multiphysics simulation methods include tightly coupled methods, which achieve high fidelity by solving the entire system jointly (Knoll & Keyes| 2004 Yu et al. 2025). often incur pro- hibitive complexity and computational costs in real-world applications (Guo et al. 2025). Loosely coupled methods, which solve each physics field separately and iteratively exchange information un- til convergence, are more widely adopted, balancing practicality, efficiency, and acceptable precision (ööhner et al. 2006). However, both approaches suffer from requirements of extensive interdisci- plinary expertise for solver development, and high computation demands (Fan & Wang 2024).
Surrogate models (MeCabe al. 2024) and neural operators (Li et al., 2025). therefore, have emerged as an alternative to accelerate multiphysics simulation, given the limitations of numeri- cal methods. Most of these models, however, rely heavily on coupled solutions as training data,
"Equal contribution. Work done intern Westlake University. 'Corresponding author.

---

Published as a conference paper at ICLR 2026
making data acquisition at least five times more computationally expensive compared to using de- coupled data (Degroote et al. 2008 Causin et al. (2055) To enable training on more accessible data from decoupled physics, surrogate models further borrow ideas from numerical solvers, such as the
Gauss-Seidel framework (Milaszewicz) 1987 and perform ADMM-like (Deng et al. 2017) iterative inference to approximate the coupled solution (Lyu et al. 2025 Gao & Jaiman (024). Although improving efficiency and leveraging decoupled data, surrogate models struggle with complex spa- tiotemporal dynamics due their limited ability to capture high-frequency, high-dimensional. and stochastic behaviors (Liang et al. 2024). The defect motivates emerging efforts in generative sim- ulation (Fotiadis et al. (025). Notably, existing generative approaches primarily focus on single- physics problems (Sun et al. (023) or directly learn multiphysics solutions from coupled data (Park) al. (02T). largely overlooking the challenge of learning coupled physies from decoupled training data. A recent study (Zhang et al., 2025) explored embedding coupling iterations into every de- noising step of diffusion models to enable coupled sampling. However, the study lacks a rigorous theoretical foundation guarantees of reliability. Emerging generative approaches show promise for high-fidelity modeling but typically rely on coupled data or lack theoretical guarantees of relia- bility. The above challenges raise critical research question: How that learns coupled physics from decoupled training data while ensuring highfidelity efficiency, and reliability?
develop a framework
To address the challenges, we introduce Generative Coupled Physics Simulation (GenCP), a novel and principled paradigm for multiphysies simulation, which achieves decoupled training and cou- pled inference with "3H", combining high fidelity from generative modeling, high efficiency from our "coupling in flow" design, and high reliability guaranteed by numerical theory.
We begin by clarifying that the essence of solving coupled problems can be reformulated as mod- eling probability density evolution in functional spaces of physical fields. During training on de- coupled data, we employ flow matching (Lipman et al., 2022) to learn the evolution of conditional probability, where each physical field evolves from noise to its target conditioned on the other one. During inference, we apply operator splitting in the functional space at each flow step, thereby merging conditional probability transitions into a coherent joint distribution transition. Our method physically corresponds to iteratively solving the coupled fields in the noisy latent space they evolve toward the solution. Figure I illustrates the overall schematic of the GenCP paradigm.
Starting from the continuity equation of the joint distribution, we show that the flow matching (Lip- man et al. 2022) sampling process can be reinterpreted with an operator operator-splitinng scheme for proba- bility density functions. We prove in Hilbert space that our proposed joint inference method enjoys error controllability guarantees. On top of the theoretical foundation, we examine GenCP on both naive and complex settings. First, we design synthetic dataset in 21) space demonstrate the joint distribution sampling effect from learning the conditional distribution. Then, we evaluate GenCP on three challenging multiphysics settings and compare against surrogate-based and other generative paradigms. Across these tasks, GenCP consistently outperforms existing approaches with error re- ducing ranging from 12.54% 42. 42.85% and even exceeding 65% on certain metrics, demonstrating its effectiveness and robustness in coupled scenarios.
The contributions of our work include: (1) We reformalize generative coupled simulation by mod- eling it as probability density evolution in the functional space of the physical field, offering prin- cipled view of the problem.
(2) We propose GenCP, an efficient generative paradigm that achieve coupled inference with decoupled learning by introducing operator splitting mechanism into the flow matching process. (3) We establish a theoretical guarantee by connecting the flow-matching continuity equation with operator splitting, and prove that our scheme enjoys controllable error bounds. (4) We conduct concrete evaluations on both synthetic setting and challenging multiphysi- cal scenarios, where GenCP consistently outperforms surrogate and generative baselines.
RELATED WORK
In this section, we briefly review the related work on both numerical and learning-based coupling simulation, highlighting key insights as well as the gaps that motivate our approach. A full discus- sion is provided in Appendix A
2

---

Published as a conference paper at ICLR 2026
Training with data of decoupled physical fields
Inference through iterative updating each physical field using the denoised other physica! fields.
Stnadare Ig beclborderr
Denoising step for fluid field
condioe
fo~N (1,
Fowr Nawhing (1180
Phytical timect
fo
pf
for
fi-at
pf
f
Tluid reld as foad boundarrs conditgn for aructtrn
90 (4, 02)
go
Ubc
1
(G.
Pmyiicl
Dince
Denoising step for displacement field of structure
Figure I:
The schematic of GenCP. With a pretrained model with decoupled physics, GenCP could manage to achieve coupled physics during flow steps. Note that here we only use Lie-Trotter Split- ting demonstration.
Numerical simulation. Classical numerical solvers typically categorized into tightly coupled (monolithic) and loosely coupled (partitioned) schemes. Monolithic approaches enforce interface conditions exactly and ensure strong stability, but are computationally prohibitive (Ruan al., 2021] Sun et al. 2024). Partitioned approaches are more efficient and easier to implement, though they are often prone to unstable convergence (MacNamara Strang 2016)] Yue Yuan 2011 Ha al.) 2023 Fourey et al.] (2017). To mitigate this trade-off, techniques such as operator splitting (Go- dunov& Bohachevsky 1959 Strang 1968; Trotter 1959) offer theoretically grounded compromise between strictly monolithic and partitioned approaches.
Learning-based simulation. When coupled training data are available, neural operators have been extended to directly model multiphysics PDEs with enhanced cross-field communication, underscor- ing the importance of functional priors (Rahman et al.) 2024 Li et al., 2025 Kerrigan et al., 2023). In settings where only decoupled data are accessible, surrogate models aim to reconstruct coupled solutions via iterative inference, yet often fail to capture high-frequency, high-dimensional dynam- ies (Sobes al. 202T; Tiba et al. (025). More recently, generative approaches have been explored to improve fidelity, but existing efforts either remain confined to single-physie scenarios lack rigorous thcoretical guarantees, limiting their applicability (Zhang et al. 2025 Sun et al. 2023).
3 METHOD
We organize the methodology into three parts. Section 3.1 formulates the coupled problem, intro- duces the weak continutty-equation perspective, and explams how decomposing the velocity field naturally leads to operator splitting as the basis for our method. Section 3.22 describes the time- parameterized linear interpolation used produce instantaneous-vellciiyy targets and the resulting training objectives. Section 3.3 presents coupled inference via operator splitting together with an informal error bound. Finally, Appendix[B] Oppendix/Bprovidde provides the precise assumptions and complete technical proofs of stability, consistency, and convergence.
3.1 PROBLEM REFORMULATION
Let D C Rd be bounded open domain where the physical fields live. The first field is a function D Rdr, F, and the second field is D - Rd 2= G. For concreteness and implementation one may take F (D2 Rdy) and G 1.2 Rd) which are separable Hilbert spaces. The joint state
(f.g) in the product space 24 F x G, I/|l || ll 1g- We consider probability laws μ on 24 with finite second moment, denoted μ E (2))). The coupled evolution of the joint state described by family (eeheeo...) In flow-matching, we seek a time- dependent velocity vector field v(t, u): [0. 1] U U, such that {ur} is transported by v.
Fully coupled trajectories are typically unavailable. Instead we assume access to decoupled single- field solver data over unit physical time step, where subseripts 0 and 1 denote the beginning and
3

---

Published as a conference paper at ICLR 2026
end of that solver step in physical time (not to be confused with the normalized interpolation time
t E [0, 1] used later for flow matching): Dj = {(fu. 90) + (fi.0)).2 D3 = {(fo. go) + (f 91)} Here 9 means that the g-field held fixed while f evolves from fo to f1. and f means that the f-field is frozen while evolves from 90 to 91. The two-field presentation is for clarity; GenCP can extend to m fields by learning from decoupled datasets where each field evolves while others are frozen, and composing them in eyclic order.
Weak continuity equation. We first recall the strong-form continuity equation. If the joint proba- bility law 14 admits a density Pe (u) on the functional state space U F x G. and if Pt is sufficiently regular, then its evolution under a time- denenden velocity field v(t, u) satisfies the strong- form-formmmo PDE
depc(u) + V.. (pe(u) v(t, u)) 0,
(I)
expressing conservation of probability mass it is transported by v. In finite dimensions this PDE is well-defined, but for distributions over function spaces several obstacles arise: empirical measures coming from decoupled physics solvers are typically singular and do not admit density P;: the divergence operator V. (prv) is not meaningful when u represents functions rather than vectors; and even when densities exist, numerical or learned models cannot reliably approximate derivatives in infiniteddimensione spaces.
These limitations motivate the of the weak continuity equation Kerrigan et al. 2023 which defines the evolution of bt only through its action on smooth test functions : 22 - R:
S.' L, ((ro())((Db())))) (Di(l), vlt, w))u) duc(u) dt = 0,
where vtt, u) is the time- deeendependent velocity vector field. Instead of differentiating the measure or its density, the weak form enforces that all observables evolve consistently with the velocity field. This provides mathematically well-posed description of measure transport in infnite-dimensionn settings such as F x G.
Decomposing the weak equation. To exploit the decoupled datasets we decompose Because the weak formulation is linear in we obtain
S. (0ip(u) + (2)0))&)(0))()))) due(u) dt = 0,
where the component Liouville operators are
(t())(( (ul u), Diplu)u
(t(t)((u) = (olo) (). Dipu))u-
Thus if we can learn By and "g from D, and Dg, their sum recovers the generator of the joint weak evolution. The decomposition enables us to consider the weak continu- ity equations formed separately for each component velocity field. For the f-component: So Su (Drp(u) + (Di(u), vlf) u))u) (u)(u) dt = O. For the g-component: So Su (Drp(u) + (Di(u), olo) w))s)) dur(u) dt = 0.
Lie-Trotter splitting. The Lie-Trotter splitting scheme provides a principled mechanism for re- combining the two learned partial dynamics during inference. Because e o(/) + v(s), the full coupled evolution can be approximated over small time step T by first evolving f with g held fixed (the flow induced by o()),, and then evolving with f held fixed (the flow induced by u(s)). In the limit of small т, this alternating update is consistent with the joint flow generated by v(f) + ((0), which explains how separately learned conditional velocities can be merged into a coherent coupled inference procedure.
To formalize this, let S.pre= denote the observable propagator (S)+he1P)(u) == ((Petter(u)) where he+hee is the flow map of the full velocity field v. For sufficiently small h. the propagator admits the first-order approximation
S
Stheet s(i) St-
((9)
4

---

Published as a conference paper at ICLR 2026
where 5() and S(9) correspond to the flows induced by vlf) and vاا), respectively. By duality, this translates to the measure-leve operator splitting
Me+h
(Ø (
(4hh-1) # Mec.
(3)
Iterating this update with step size yields the classical Lie-Trotter product formula used during inference (Trotter 1959] Strang 1968 Bätkai et al. 2011). This justifies recombining the separately learned conditional velocity fields through alternating partial flows, thereby recovering an approx- imation of the true coupled dynamics. Rigorous stability, consistency, and convergence results are deferred to Appendix
3.2
TIME-PARAMTTERIEED LINEAR INTERPOLATION AND TRAINING
We now give fully constructive procedure that produces instantanecous-veloiii targets from de- coupled data.
Reference distributions. Choose tractable reference measures T& F and Ti on G. In practice, these may be simple Gaussian priors defined on the finite discretization of D.
Learning the velocity for the field. Sample (fi,g) ~ Dr and independent references 2, T5, 2g ~ Tg. For & [0, 1] define the linear interpolants ft = (i - t)zf + tf1 and 9t = (1 -)2g + tg.
The instantaneous derivative along this interpolation vector
dfc f di
We use (fe, St. t) as input to the conditional velocity model for ff and take as target. This setup yields random supervisory signal whose expectation conforms the conditional velocity stipulated in the weak continuity equation. Symmetrically, sample (J (1) ~ D, with independent references z'r1 dgg and set ft (1 - t) + tf. = (1-t)z' + tgg. The vector is given by Then (fr, 9v, is the input and by is the target for the conditional velocity of g.
dg: dt Ug 91 2''
Learning objectives. Parameterize two operator-valuce models êy t; 0y) and (f g, t; 0g) that map inputs (fi, Dt, 1) to elements of F and 9, respectively. Minimize the -quare losses
()(0,) B),((r [|ws - ff 9t. t; t; 11
(4)
(0y0 ).((f=)))Dp)) [lle - fy 9t, t; 1.)]
(5)
The two training procedures use different decoupled datasets, which match typical solver outputs. The specific training algorithm can be found in Appendix F.
3.3 COUPLED INFERENCE VIA OPERATOR SPLITTING
Coupled inference. After training we obtain learned conditional velocity estimators úy and ég- We then recombine their induced flows using Lic-Trotter splitting to produce coupled samples. Concretely, let A be the flow map obtained by integrating ((4) = (ūy, 0) and ( the flow map for (B) = (0, 4²). For step size 1/N, the ieeTTroteer push-forward update is given by Eq. which at the state level corresponds to the explicit alternating update:
Algorithm 1 Coupled Sampling via Lie-Trotter Splitting
Input: initial uo (fo. 9o). trained operators if, 0g. steps N
2: u t
3: for = 0 to N
do t k/N
4:
5:
f f + à
6:
g (f
7: u (f
8: end for
9: return
5

---

Published as a conference paper at ICLR 2026
Theoretical guarantee of acceptable error bound. The two sources of the learned split- ting scheme are the splitting discretization and the learning approximation. Under the stability and regularity hypotheses stated precisely in Appendix B one can show the following bound.
Theorem 3.1 (Informal erorr bound), Let μ1 denote the true terminal measure at l t= 1 generated
by the unknown joins velocity v. Let pi (r.eam)) denote the terminal measure obtained by applying the learned Lie-Trotter composition with step size using 0g 0g. If the learned fields uniformly approximate the true conditional velocities on the region visited by the dynamics with errors at most et Eg. then
Wi (1( (T.learm) .ui) S Cstab (T + 5 +
where W1 is the Wasserstein-1 distance and Cstab depends on Lipschitz and growth constants for the flows.
The term J comes from the first-order nature of Lie-Trotter splitting, meaning the error introduced by time-splitting scales linearly with the step size. The terms Eg,Eg capture how well the learned conditional velocities approximate the true conditionals. Thus, improving either time resolution or regression accuracy reduces the total error.
Precise assumptions, the stability lemma for propagator composition, the infinitesimal consistency of the split generator, and the detailed convergence proof are collected in Appendix B Those techni- cal statements track which Lipschitz and growth conditions are needed and make the informal bound above rigorous.
4 EXPERIMENTS
In this section, we aim to answer the following 3 questions: (1) In the purely mathematical setting of probability distributions, can GenCP enable training two separate flow matching models on datasets describing conditional distributions, and then directly combine them at inference sample from the joint distribution? (2) Beyond the mathematical case, can GenCP maintain strong performance (effectively modeling coupling behaviors at inference with only decoupled training) in complex and high- high-dimensional coupled simulation problems? (3) In addition to simulation accuracy, how does the inference efficiency of GenCP compare with existing paradigms? Next, question (1) is addressed
in Section 4.2 and questions (2) and (3) are answered in Section 4.3 and Section 4.4
4.1 EXPERIMENTAL SETUP
We evaluate GenCP on four representative cases: a simple 2D distribution, designed to directly illus- trate our conditional-to-joine sampling paradigm, two canonical FSI benchmarks (Double-Cyline) and Ture-Hron) and a complex nuclear-thermal setting. To assess the capability of GenCP, we compare against two alternative paradigms for modeling coupled physics, each with two neural operator backbones (CNO and #N The two alternative paradigms represent (i) surrogate model-based Picard-iteration approach and (ii) the M2PDE paradigm.
Across all tasks, training is conducted exclusively on decoupled datasets, following Algorithm In practice, this means we model the distributional evolution of one variable while treating the other as a known control condition. A subset of the decoupled data is reserved as a validation set to evaluate learning on the decoupled setting, and the coupled data serve as the test set to assess the coupling performance. Further details on the datasets used can be found in Appendix [
At inference time, GenCP performs coupled simulation by alternately evolving from noise to the target distribution using the learned conditional distributions via operator splitting. Specifically, we adopt the Lie-Trotter splitting method, as shown in Algorithm T to achieve a balanced trade-off between precision and computational cost.
For comparison among the paradigms, the surrogate- based paradigm trains a neural network that end-to-end maps system inputs to outputs decoupled data and then applies Picard iteration at inference to approximate the coupled solution. The M2PDE paradigm instead trains a diffusion
ONNO* denotes the FNO framework with SiT as its internal model, making the backbone more expressive for learning objectives while retaining its functional properties.

---

Published as a conference paper at ICLR 2026
Table 1:
Comparison of statistical distances between the target joint distributions and those gener-
ated from different paradigms.
Easy Distribution
Complex Distribution
GenCP
M2PDE
GenCP
M2PDE
W
0.4366
0.5177
0.4928
25450.5442
MMD
0.0095
0.0141
0.0053
inf
Energy Distance
0.0411
0.0625
0.0061
332.3619
Original
M2PDE
GenCP
Samples
Figure 2: Visualization sampling results on complex distribution with different paradigms the synthetic setting. In order to examine behavior except outlier sampling, the visualized results of M2PDE also include a zoomed-in in view on the same coordinate region as the other views.
model on each decoupled dataset, using the DDIM estimation the denoised result of one field as the condition to denoise the other. Detailed implementation settings and pseudo-code for baseline paradigms are provided in Appendix
Full implementation details (including backbone architectures, sampling paradigms, and training configurations) are provided Appendix. As an additional contribution, we also open-source the dataset used in this study to support further method development in coupled simulation.
4.2 SYNTHETIC SETTING ON 2D DISTRIBUTION
To intuitively demonstrate the core idea of GenCP, that learning p(>: y) and p(y x)) on decoupled samples suffices for direct joint sampling of p(x,, y) during inference, we first investigate it on the synthetic dataset. We note that M2PDE (Zhang et al. 2025) also tried to sample on the joint dis- tribution with generative model based on training with conditional data. So here compare our performance with M2PDE.
We design two synthetic distributions synthetic dataset: an "Easy distribution" based simple Gaussian mixture, and a "Complex distribution" with multiple Gaussian components arranged in multimodal pattern. To obtain conditional information while avoiding direct joint sampling, generate two complementary datasets for training that separately focus on p(y ( and p(zly) with perturbed marginal for x and y, respectively. This setup ensures that models are trained only on conditionals, while the underlying target joint structure remains nontrivial. Further details of the synthetic setting could be seen in Appendix D
We visualize the sampling effect in Figure and report quantitative metrics, including Wasserstein-1 distance, MMD, and Energy distance between the target distribution and the generated one Table I The distributional distance results show that our paradigm slightly outperforms the bascline on simple distribution and achieves substantial advantage on more complex distribution, maintaining stable performance while the baseline degrades severely. These results demonstrate that GenCP can successfully recover the true joint distribution with only conditional training, closely matching the ground truth both visually and in distributional metrics. In contrast, the baseline method performs poorly: not only fails to capture the main body of the distribution, but also exhibits strong instabil- ity that leads to frequent sampling of outliers. This weakness arises from its design, as well as the error accumulation induced from using intermediate estimates as conditions, ultimately causing mode collapse drift, especially in complex distributions.

---

Published as a conference paper at ICLR 2026
4.3 APPLICATION TO FSI SCENARIOS
We next evaluate GenCP two fluid-structure coupling tasks, which exemplify real multiphysics coupling with strong bidirectional feedback. Both settings are widely used benchmarks for FSI and present significant challenges: strongly coupled fields, high-dimensional and high-frequency features, as well as nonlinear feedback loops. The FSI tasks here are formulated as follows: given known steps of both fields, the goal is to predict the couple dynamics over the next 12 steps, using a model trained only with decoupled data.
In both settings, we describe fluid behavior using velocity fields (u, v) and a pressure field (p), while structural dynamics are represented by signed distance field (SDF). Thus, the multiphysics representation of the FSI problem can be composed into a four-channel field defined on the same computational domain. Through this careful design, we avoid the common difficulty encountered when representing fluid-structure interface in normal FSI algorithms, such as Arbitrary Lagrangian- Eulerian (ALE) (Korobenko et al. 2018) and Immersed Boundary Method (IBM) (Tian 2014).
Turek-Hron. Turek-Hron is a classical FSI case: rigid cylinder is placed in cross-flow, followed by flexible beam that bends under the influence of vortices shed in the wake. This setup is also widely used as a standard for evaluating numerical FSI solvers, making it a stringent test of fidelity. We conduct training and inference following the procedure described the previous section.
Table 2 compares the relative L2 norm errors and inference costs of the three inference paradigms across two backbone models. On both FNO* and CNO backbones, GenCP consistently achieves bet- ter test accuracy on coupled data than both surrogate- baarroggte-baeeddddeee and M2PDE paradigms. With the FNO* backbone, the average error across four fields about 26.77% lower than the best-performing base- line, and with the CNO backbone, the error is about 12.54% lower. Furthermore, the CNO backbone even achieves results on SDF that are close to those of Joint Training. paradigm demonstrates extremely pronounced advantage: thanks to the operator-splitiine design that embeds coupling directly into the flow matching process, GenCP requires only 10 sampling steps to generate accurate coupled solutions. This design allows GenCP to simultancously maintain accuracy while significantly improving inference efficiency compared to surrogate- based explicit iterations diffusion-based iterative coupling condition estimated fields.
terms of efficiency, our
Table 2: The relative 1.2 norm error and inference time of different methods the Turek-Hrow setting.
L.2 Norm Field
Validation on Decoupled Data
Test on Coupled Data
Inference Time
SDF
14
SDF
Joint Truining
0.0088
0.03441
0.0544
0.0079
MOPDE-FNO*
0.0452
0.1582
0,1444
0.0206
0.0590
0,2415
0,2474
0.2482
277.20s
Surrogate-FNe FNO*
0.0181
0.0756
0.0953
0.0087
0.0550
0.2257
0.2553
0.0112
93.20s
Our GenCP-BNO*
0.0091
0.0471
0.0548
0.0069
0.0396
0.1678
0.1897
0.0081
19.50s
M2PDE-CNO
0.0497
0,2075
0.3466
0.039377
0.05769
0.2809
0.4087
0.0390
347.00s
Surrogute-CNO 0.0204
0.0994
0,1526
0.0205
0.0469
0.1888
0.2278
0.0242
300.25s
Our GenCP-CNO
0.0150
0,0626
0.1187
0.0152
0.0388
0.1821
0.2166
0.0183
16.25s
Figure 3 provides visual comparison of the coupled solutions for each field across different paradigms using FNO* as the backbone. As shown, our method closely matches the ground truth, with error maps confirming its ability to preserve both global patterns and the high-frequency details induced by fluid-structure interactions. and it clearly outperforms the bascline methods. In partic- ular, examining the reconstructed SDF field highlights that, although the surrogate- method reports relatively low error values in Table it fundamentally fails to model the oscillatory bending dynamics of the beam. By contrast, our GenCe,while trained only on decoupled data, is the only method can can successfully capptuu these bending sendinggefffects. which are intrinsic to true coupling. In fact, the seemingly low error of surrogate models arises from their deterministic end-to-end predic- tion. In contrast, our probabilistic modeling approach accounts for both mode errors and stochastic noise, explaining why our visual results appear far superior to the basclines, even though the quan- titative metrics show modest improvements.
8

---

Published as a conference paper at ICLR 2026
Surrogate- ENO
MEPDE-FNO
GenCP-ENO
Figure 3: The visualization of results from different paradigms using FNO* as backbone on Turek- Hron setting. (Note that the structural boundary shown in the fluid field visualization is masked from the ground truth and does not represent actual structural deformation.)
Table 3: The relative L2 norm error of different methods on the Double Cylinder setting.
Rel L2 Norm Field
Validation on Decoupled Data
Test on Coupled Data
SDF
SDF
Joint Training
0.0110
0.0387
0.0786
0.0034
M2PDE-FNO* 0,0128
0.0414
0,2238
0.0065
0.0571
0.2583
0.8879
0.0152
Surrogate-NNO=
0.0189
0.0625
0.0681
0.0041
0.0717
0.3058
0.4979
0.0196
Our GenCP-FNO*
0.0098
0.0441
0.0525
0.0060
0.0522
0.2397
0.3987
0.0061
M2PDE-CNO
0.0295
0.2072
0.6468
0.0156
0.0379
0.2694
0.9966
0.0179
Surrogate-CNO
0.0364
0.1078
0.2060
0.0067
0.0612
0.2058
1.0722
0.0080
Our GerCP-CNO
0.0141
0.0474
0.1753
0.0058
0.0279
0.1150
0.6208
0.0055
Double cylinder. To evaluate the performance of our method on more strongly coupled problems, we construct more challenging FSI scenario by increasing the time step. Specifically, we adopt a double-cylinde configuration, in which one cylinder is fixed in a cross-flow while the other is attached to a linear spring and allowed to oscillate vertically with respect to inlet. This setup in- duces strong bidirectional feedback between vortex shedding and structural oscillation. The training and inference procedures are identical to those described previously and are not repeated here.
The error analysis of GenCP compared with bascline paradigms is quantified in Table 3 The trends are consistent with those observed the Turek-Hron case. Evidently, our method significantly out- performs baseline methods across all physical fields, with average 34.44% on FNO* backbone and up to 42.85% on CNO backbone. These results demonstrate that GenCP exhibits
reductions of approximately
greater advantages over baselines in strongly coupled problems.
In addition, we also report results obtained by training flow matching models directly on coupled data, iee., learning the joint distribution. Theoretically, this approach represents the upper bound

---

Published as a conference paper at ICLR 2026
of our method. In other words, the discrepancy between the results of joint training and those of training on decoupled data reflects the error introduced by "conditional learning to joint sampling". Although the error values from joint training remain lower, the performance of GenCP is already very close to this upper bound compared to other baselines, highlighting the value of our proposed paradigm. For additional visualization results and analyses, please refer to Appendix E
4.4 APPLICATION TO NT COUPLING SCENARIOS
To further demonstrate the scalability and generality of GenCP, we include an additional experiment on a Buclear-Thermal Coupling problem. Here GenCP can be seamlessly applied to this more complex multiphysics scenario, involving scaling to additional physical fields, and stronger or more intricate coupling mechanisms.
The NT coupling system requires jointly solving neutron diffusion, solid heat conduction, and fluid heat transfer equations, while accounting for (i) the negative thermal feedback between neutron physics and temperature, (ii) the unidirectional influence of the fluid field the neutron field through temperature, and (iii) the strong interface coupling between the fuel and the coolant. The geometry and coupling configurations of this system are illustrated in Figure [ In this setting, our goal is to predict the system's transient evolution of coupling physics under different boundary conditions (neutron flux and fuel temperature measurement at the left side), with learning only on completely decoupled data.
The visualizations of the experimental results with three paradigms are provided in Figure As shown in the Table GenCP consistently outperforms both baselines across the neutron, fluid, and solid fields. Specifically, using FNO* as the backbone, GenCP achieves an average error reduction of 49.9% relative to M2PDE and 51.7% relative to the Surrogate paradigm. With the CNO backbone, the improvements are even more substantial, with error reductions of 58.2% over M2PDE and 78.8% over Surrogates. These results demonstrate that our method significantly surpasses existing baselines tasks in complex multiphysics settings.
on "decoupled-training to
To further investigate the source of errors, we also report the validation error on the coupled dataset, obtained by models evaluated using the ground-truth auxiliary coupled fields as known conditions. This eиror corresponds to the approximation error in Theorem 3.1 while the discrepancy between this prediction and the actual coupled inference result reflects the splitting error т. The results indicate that, consistent with theoretical guarantees, the error introduced by our conditional-to- joint sampling method is minimal.
Table 4: The relative 1.2 norm error of different methods on the Nuclear Therraal Coupling setting.
Rei 12 Nor
Decoupted validation on decoupled data
Decoupled validation on coupled data
Coupled test coupled data
Fiehd
Neation
Fud
Fliid
Nautun
Fuid
Neut
Fluid
Our GANCP-FNO
0.0022
0.0006
0.0038
0.0081
0.0571
B.0364
0,0085
0,0364
0.1270
END*
0.0056
0,0014
0,0022
0.0140
0.0167
0.0767
0.0149
0,0576
0.1095
M2PDE-PNO
0.0022
0.0014
0.0018
0.0082
0.0085
0.0237
0.0136
0,1257
0.0463
Our GeICP-CNO
0.0024
0,0005
0,0110
0.0047
0.0083
0,0308
0,0044
0.0105
0.1330
Sumpps-CNO
0.0046
0.0007
0.0082
0.0073
0.0044
B.0667
0.0138
0.0553
0. 30306
MIP06-CNO
0.0397
5 CONCLUSION
In this work, we tackle the fundamental challenges of coupled simulation, where efficiency, fidelity, and the ability to leverage decoupled data remain major bottlenecks. To address these challenges, we reformulated the problem from functional probabilistic evolution perspective, recognizing that the essence of coupling is sampling from the joint distribution via learned conditional distributions. Building on this insight, we have introduced GenCP, principled generative paradigm that enables decoupled training and coupled inference in flow matching, achieving substantially higher accuracy and efficiency than existing methods across three distinct scenarios. Moreover, by applying operator- probabilityydde density evolution process, we established provable error bound for this conditional-ta -o-joot sampling scheme, providing GenCP with a solid theoretical foundation.
splitting theory to the probability
Looking ahead, we believe this work paves the way toward probabilistically principled and prac- tically scalable generative paradigms for coupled simulation, uniting theoretical guarantees with real-world applicability.
10

---

Published as a conference paper at ICLR 2026
ETHICS STATEMENT
This work proposes generative modeling paradigm for coupled simulation using synthetic and simulated data and provides theoretical guarantee for it. All datasets are collected without involving sensitive information.
REPRODUCIBILITY STATEMENT
The code is available at github. gethub.com//Id4e.. com/AI4Science. WeeeseeeeeeetlaakU/GGGCCPCP We provide a unified and mod- ular code framework, together with seripts for reproducing all experiments. We will also release the data and checkpoints with full documentation to ensure transparency and reproducibility.
REFERENCES
Sergio Blanes, Fernando Casas, and Ander Murua. Splitting methods for differential equations. Acta Numerica, 33:1-161, 2024.
Martina Bukač, Guosheng Fu, Anyastassia Seboldt, and Catalin Trenchea. Time-adaptive partitioned method for fluid-structure interaction problems with thick structures. Journal of Computational Physics, 473:1 11708, 2023.
Martina Bukač, Ivan Yotov, and Paolo Zunino. An operator splitting approach for the interaction be- tween fluid and a multilayered poroelastic structure. Numerical Methods for Partial Difserential Equations, 31(4), 11 2014.
András Bátkai, Petra Csomós, Bälint Farkas, and Gregor Nickel. Operator splitting for non- autonomous evolution equations. Journal of Functional Analysis, 260(7)21133-2190 2011. ISSN 0022-1236.
Sunčica Canić, Marija Galić, and Boris Muha. Analysis of a 3d nonlinear, moving boundary problem describing (9):6621-6681, 2020.
interaction. Transactions of the American mathematical society, 373
Paola Causin, Yean-Frédéric Gerbeau, and Fabio Nobile. Added-mass effect in the design of par- titioned algorithms for fluid-structus problems. Computer methods in applied mechanics and engineering, 194(42-44):4506- 4 2005.
Joris Degroote, Peter Bruggeman, Robby Haelterman, and Jan Vierendeels. Stability of coupling technique for partitioned solvers in fsi applications. Computers & Structures, 86(23-24):2224 2234, 2008
Wei Deng, Ming-Jun Lai, Zhimin Peng, and Wotao Yin. Parallel multi-block admm with o (1/k) convergence. Journal of Scientific Computing, 71 2017.
Xiantao Fan and Jian-Xun Wang, Differentiable hybrid neural modeling for fluid-structure interac- tion. Journal of Computational Physics, 496:1 12584, 2024.
Carlos A Felippa, Kwang- Chun Park, and Charbel Farhat. Partitioned analysis of coupled mechan- ical systems. Computer methods in applied mechanics and engineering, (90024-25)33247-32770 2001.
Stathi Fotiadis, Noah D Brenowitz, Tomas Geffner, Yair Cohen, Michael Pritchard, Arash Vahdat, and Morteza Mardani. Adaptive flow matching for resolving small-scale physics. In Forty-second International Conference on Machine Learning, 2025.
G Fourey, Corentin Hermange, David Le Touzé, and Guillaume Oger. An efficient fsi coupling strat- egy between smoothed particle hydrodynamics and finite element methods. Computer Physics Communications. 217:66 6--8,, 2017.
Rui Gao and Rajeev K Jaiman. Predicting quid-structurr interaction with graph neural networks. Physics of Fluids, 36(1), 2024.
11

---

Published as a conference paper at ICLR 2026
M Barzegar Gerdroodbary and Sajad Salavatidezfouli. A predictive surrogate model based on linear and nonlinear solution manifold reduction in cardiovascular fsi: A comparative study. Computers in Biology and Medicine, 189: 189:109959, 2025.
Sergei K Godunov and Ihor Bohachevsky. Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics. Matematičeskij sbormik, 47(3):271- 306, 1959.
Shihan Guo, Yizhong Sun, Yifan Wang, Xiaohe Yue, and Haibiao Zheng. fully parallelizable loosely coupled scheme for fluid-poroelastic structure interaction problems. SIAM Journal on Scientific Computing, 47(4): 2025.
Sang Truong Ha, Hyoung Gwon Choi, Ngo Cu Long, and Sang-Wor Wook Lee. A semi-implicit fi- nite element formulation of the partitioned method for ftuid-structure interaction based on flux boundary condition of pressure equation. Journal of Mechanical Science and Technology. 37(5):
2417-2424, 2023.
Ernst Hairer, Gerhard Wanner, and Christian Lubich. Oscillatory differential equations with vary- Algorithms for
ing high frequencies. In Geometric Numerical Integration: Ordinary Differential Equations, pp. 531-565. Springer, 2006.
Gene Hou, Jin Wang, and Anita Layton. Numerical methods for fluid-structure interaction a view. Communications Computational Physics, 12(2):337-377, 2012.
Jiawei Hu, Zihao Dou, and Weiwei Zhang. Fast fluid-structure interaction simulation method based deep learning flow field modeling. Physics of Fluids, 36(4), 2024.
Haochen Huang, Daogang Lu, Yu Liu, Danting Sui, Fei Xie, and Hao Ding. Research on a rapid prediction method for multi-physics coupled fields in small lead-cooled fast reactors based on machine learning. Nuclear Engineering and Design, 438: 14065, 2025.
Casey Icenhour, Shane Keniley, Corey DeChant, Cody Permann, Alex Lindsay, Richard Martineau, Davide Curreli, and Steven Shannon. Multi-physics object oriented simulation environment (moose). Idaho National Laboratory (INL), Idaho Falls, ID (United States), 11 2018. URL
https: /tt::://wwwwwsro ostii.ssti... wtt.ggiblii/ol190000
Gavin Kerrigan, Giosue Migliorini, and Padhraic Smyth. Functional flow matching, arXiv preprint arXiv:2305.17209, 2023.
Dana A Knoll and David E Keyes. Jacobian-free newton -rwtoon-krylo methods: a survey of approaches and applications. Journal of Computational Physics, 193(2):357-397, 2004.
Kazuma Kobayashi, Jaewan Park, Qibang Liu, Seid Koric, Diab Abueidda, and Syed Bahauddin Alam. When network architecture meets physics: Deep operator learning for coupled multi- physics. arXiv preprint arXiv:2507.0366h 2025.
Artem Korobenko, Yuri Bazilevs, Kenji Takizawa, and Tayfun E Tezduyar. Recent advances in ale-vms and st-vms computational aerodynamic and fsi analysis of wind turbines. Frontiers in Computational Fluid-Strucce Structure Interaction and Flow Simulation: Research from Lead Investiga- tors under F018 pp. 253-336, 2018
Shibo Li, Tao Wang, Yifei Sun, and Hewei Tang. Multi-physics simulations via coupled fourier neural operator. arXiv preprint arXiv:2501. 17296, 2025.
Aoming Liang, Zhaoyang Mu, Pengxiao Lin, Cong Wang, Mingming Ge, Ling Shao, Dixia Fan, and Hao Tang. m²m:: Learning controllable multi of experts and multi-scale operators are the partial differential equations need. arXiv preprint 1rXi::2410011617. 2024.
Yaron Lipman, Ricky TQ Chen, Heli Ben- Gen-Hamu, Maximilian Nickel, and Matt Lc. Flow matching for generative modeling. arXiv preprint arXiv:2210.022 227.02777 2022.
Rainald Löhner, Juan R Cebral, Chi Yang, Joseph D Baum, Eric L Mestreau, and Orlando Soto. Extending the range and applicability of the loose coupling approach for fsi simulations. In Interaction: Modelling, Simulation, Optimisation, pp. 100. Springer, 2006.
12

---

Published as a conference paper at ICLR 2026
Yanfang Lyu, Yunyang Zhang, Zhiqiang Gong, Xiao Kang, Wen Yao, and Yongmao Pei. A novel hybrid neural network of ftuid-structure interaction prediction for two cylinders in tandem ar- rangement. Engineering Applications of Computational Fluid Mechanics, 19(1):2513663 2025.
Shev MacNamara and Gilbert Strang. Operator splitting. In Roland Glowinski, Stanley J. Osher, and Wotao Yin (eds.), Splitting Methods in Communication, Imaging, Science, and Engineering, Scientific Computation. Springer, Cham, 2016.
ZM Malikov, ME Madaliev, SL Chernyshev, and AA Ionov. Validation of a two-fluid turbulence model in comsol multiphysics for the problem of flow around aerodynamic profiles. Scientific Reports, 14(1):2306, 2024.
Michael McCabe, Bruno Régaldo-Saint Blancard, Liam Parker, Ruben Ohana, Miles Cranmer, Al- berto Bietti, Michael Eickenberg, Siavash Golkar, Geraud Krawezik, Francois Lanusse, et al. Multiple physics pretraining for spatiotemporal surrogate models. Advances in Neural Informa- fion Processing Systems, 37:1 1970111 119335, 2024.
S Meduri, M Cremonesi, U Perego, O Bettinotti, A Kurkchubasche, and V Oancea. A partitioned fully explicit lagrangian finite element method for highly nonlinear fluid-structure interaction problems. International Journal for Numerical Methods in Engineering, 1 113(1):43-64, 2018.
Juan Pedro Milaszewicz. Improving jacobi and gauss-seidel iterations. Linear Algebra and Its
Applications, 93:161-170, 1987.
Tae-Ryong Park, Hyunseong Park, Kiyoul Kim, Chae-Nam Im, and Jang-Hyeon Cho. Heat and weight optimization methodology of thermal batteries by using deep learning method with multi- physics simulation. Energy Conversion and Management, 236:1 14033, 2021.
RWCP Pramaník, RWCP Verstappen, and PR Onck. Computational fluid strucutre interaction in biology and soft robots: A review. Physics of Fluids, 36(10), 2024.
Md Ashiqur Rahman, Robert Joseph George, Mogab Elleithy, Daniel Leibovici, Zongyi Li, Boris Bonev, Colin White, Julius Berner, Raymond A Yeh, Jean Kossaifi, et al. Pretraining codomain attention neural operators for solving multiphysics pdes. Advances in Neural Information Pro- cessing Systems, 37: 37:104035-104044 2024.
Liangwang Ruan, Jinyuan Liu, Bo Zhu, Shinjiro Sueda, Bin Wang, and Baoquan Chen. Solid-fluid interaction with surface-tension-dominnn contact. ACM Transactions on Graphics (TOG), 40(4): 1-12, 2021.
Vladimir Sobes, Briana Hiscox, Emilian Popov, Rick Archibald, Cory Hauck, Ben Betzler, and Kurt Terrani. Ai-based design of nuclear reactor core. Scientific reports, 11(1):19646, 2021.
Gilbert Strang. On the construction and comparison of difference schemes. SIAM journal nu- merical analysis, 5(3):506-517 1968.
Luning Sun, Xu Han, Han Gao, Jian-Xun Wang, and Liping Liu. Unifying predictions of deter- ministic and stochastic physics in mesh-reduced space with sequential flow generative model. Advances in Neural Information Processing Systems, 36:60636-60666 2023.
Zixian Sun, Zhixin Zeng. Jiasheng Li, and Xiong Zhang. An immersed multi-matend arbitrary lagrangian-culerian finite element method for fluid-strrucuure-nttraciioo problems. Computer Methods in Applied Mechanics and Engineering, 432:117395 2024.
Fang-Bao Tian. Fsi modeling with the dsd/sst method for the fluid and finite difference method for the structure. Computational Mechanics, 54(2):581-589, 2014.
Azzeddine Tiba, Thibault Dairay, Florian De Vuyst, Iraj Mortazavi, and Juan Pedro Berro Ramirez. Machine-learmnn enhanced predictors for accelerated convergence of partitioned fluid-structure interaction simulations. Computer Physics Communications, 300:10 2025.
Hale Trotter. On the product of semi-grouph of operators. Proceedings of the American Mathe- matical Society, 10(4):545-551, 1959.
13

---

Published as a conference paper at ICLR 2026
Jianguo Wang. Kai Wei, and Haowei Cai. Dynamic response of cylindrical bridge pier subjected to the impact of flash flooding considering the effect of soil-structure interactions. In Structures, volume 72, pp. 108199. Elsevier, 2025.
Gabriel D Weymouth and Dick KP Yue. Boundary data immersion method for cartesian-grid simula- tions of fluid-body interaction problems. Journal of Computational Physics, 23:(66)::233 2011.
Kasun Wijesooriya, D Mohotti, A Amin, and K Chauhan. An uncoupled fluid structure interaction method the assessment of structural responses of tall buildings. In Structures, volume 25, pp. 448-462. Elsevier, 2020.
Changfan Yang, Lichen Bai, Yinpeng Wang, Shufei Zhang, and Zeke Xie. Multiphysics bench: Benchmarking and investigating scientific machine learning for multiphysics pdes. arXiv preprint arXiv:2505. 17575, 2025.
Y Yu, H Park, A Novak, and E Shemon. High fidelity multiphysics tightly coupled model for lead cooled fast reactor concept and application statistical calculation of hot channel factors. Nuclear Engineering and Design, 435:1 13915, 2025.
Jingyan Yue and Guangwei Yuan. Picard-newton iterative method with time step control for multimaterial non-equilibrium radiation diffusion problem. Communications in Computational Physics, 10(4):844 866, 2011.
Tao Zhang. Zhenhai Liu, Feipeng Qi, Yongjun Jiao, and Tailin Wu. M2PDE: Compositional gener- ative multiphysics and multi-component PDE simulation. In Forty-second International Confer- ence on Machine Learning, 2025.
A RELATED WORK
A.1 NUMERICAL COUPLING SIMULATION
Traditional numerical simulation for coupling physics mainly target at solving major difficulties about ensuring accuracy, stability, and efficiency (Felippa al. [2001] Hou et al., 2012). Classical numerical solvers distinguish between tightly coupled (monolithic) and loosely coupled (partitioned) schemes. Tightly coupled solvers assemble fluid and solid equations into a coupled system, enfore- ing interface conditions exactly and offering strong stability, but at the expense of high computa- tional cost. |Raan et al. (2021) introduce monolithic framework that accurately handles momentum transfer at liquid-solid-air interfaces. [Snn et al. (2024) proposed an immersed multi-material arbi- trary -aggaangian-Fullrian method to realize monolithic coupling framework.
Loosely coupled solvers, in contrast, partition system into subproblems for more practical solution, while they may suffer from unstable convergence. Besides, equipped with numerical techniques such as operator splitting (MacNamara & Strang 2016) and Picard iteration (Yue & Yuan 201T), they typically enjoy theoretical guarantees on acceptable error bounds. Representative developments include improvements in numerical integration strategies, such as semi-implicit coupling schemes (Ha et al. 2023) and time-adaptive partitioned methods (Bukač et al.) 2023), as well as advances in representation techniques, such as coupling SPH with FEM (Fourey et al. 2017), and a partitioned explicit Lagrangian FEM tailored for highly nonlinear multiphysics problems (Meduri et al.) 2018).
More specifically about operator splitting, it provides theoretically grounded compromise between strictly monolithic and partitioned approaches. By decomposing the global operator into sub- operators advanced sequentially, splitting schemes exploit sub-solvers while maintaining control over stability and eтror (Godunov & Bohachevsky 1959 Strang 1968 Trotter 1959). Extensive surveys have confirmed that Trotter and Strang compositions remain central tools in multi-
physics applications (Blanes et al. 2024 Bukač et al. 2014 Čanié et al. 2020
A.2 LEARNING-BASER SIMOLATION
Existing learning-basse based approaches for multiphysics simulation have also basically evolved along two categories: tightly coupled (monolithic) and loosely coupled (partitioned) methods. The first
14

---

Published as a conference paper at ICLR 2026
category leverages neural operators to directly simulate multiphysics PDEs. CANO Rahman al. 2024) introduced codomain attention to enhance information exchange across coupled PDE outputs, demonstrating its potential for improving gemeralizability in multiphysics problems. COMPOL (Li et al. 2025) further incorporated explicit cross-field interaction mechanisms and gained notably increased accuracy. Kobbaashi et al. (2025) systematically examined the relationship between op- erator architectures and the coupling strength of multiphysics. Collectively, these works underscore the importance of embedding functional priors (Kerrigan al. 2023) for multiphysics simulation.
The second category employs Al models to learn decoupled physics in systems and then recomposes them to obtain coupled solutions (Sobes al., (00I). In nuclear engineering, Al methods have been proposed for modeling of complex thermo- Mhermo-hydraulic hydraulic- neutronic interactions (Huang et al. [2025.. For unsteady problems, Tiba al. 2025)) incorporated reduced-order models for enabling better convergence and efficiency. These studies demonstrate the effectiveness of surrogate models in ac- celerating coupled simulations while keeping accuracy. More recently, generative models have also been explored for coupled simulation (Zhang al. 2025), like applying compositional generation to capture joint dynamics across multiple conditional models. However, its lack of theoretical guar- antees limits its applicability in strongly coupled scenarios.
B DEFINITIONS AND PROOFS
B.1 OPERATOR SPLITTING SCHEME: ASSUMPTIONS AND CONVERGENCE GUARANTEES
This appendix provides the operator-spliting analysis underlying the error control bound in Sec-
tion 3.3
To ensure the well-posedness of the dynamics and convergence of the splitting scheme, impose the following:
Assumptions.
(AI) Local Lipschitz and growth. For every R > O, there exists LR > 0 such that for all E [0, 1] and u, ū E l with [|u]ee, ||ūller S R,
||u(t, u) - utt, u)ll S LRllu- ullu.
|| w)llu S Cillulle + Cz.
This ensures local existence and uniqueness of solutions.
(A2) Well- rell-posed of subproblems. For every initial state U, the A- and B-subproblems (evolving f with g frozen, and y with f frozen) are well-posed, generating evolution fami- lies U (1, and (t( They satisfy:
(i) Existence: for any starting time 8, there exists a dense subspace Y, on which classical solutions exist.
(ii) Uniqueness: each y E Y, generates unique solution.
(iii) Continuous dependence: solutions depend continuously on (s, v). uniformly on com- pact intervals.
(iv) Exponential boundedness: there exist constants M > 1, E RR such that
[t, < Me
||lll.
V 2 s.
(A3) Accuracy of learned vector fields. On the bounded state set visited during evolution, the learned vector fields é i satisfy uniform approximation bounds
sup ||uy(t, u) éy (1, u) S E wEB
sup (t, (t, u)lla S Eg.
wEB
and are Lipschitz in u with a uniform constant T-mode)-
Lie-Trotter splitting. nuity equation
With = ((4) + BB) define the Liouville operators from the weak conti-
(a())pu))= (1)). D Dy((1))FF (B())((w) = (1, Do
so that C(t) = Ct(t) + (t). The splitting scheme evolves measures as
(1)) (Φ (B)
(
(T)
== 1/1nn.
15

---

S(A) Statt te and (B) Siuittte satisfies the Published as a conference paper at ICLR 2026
where ()) and (P)) denote flows generated by the learned fields iy 0g. The global approximation over [0, 1] is the n-fold composition
(7) ( $b(1) Mo.
Stability. Assume there exist constants corresponding propagators S(4) ((())
2 w E R such that for any partition {4x} and
III
S B) 1
I Me",
C(
for Banach space B of observables controlling sup-norms and gradients.
Consistency. The scheme is consistent:
S(4) 1 L(t)p
Sithr
lim
h
for any yp T. This follows from the Taylor expansion of the A- and B-flows on finite-dimensional projections.
Convergence.
The local splitting error
S(A) (4t p(u) L(t)p(u),
Rx u)
which Ohh) under regularity conditions. Thus, the global error first order O(7) when using exact vector fields. With learned fields, an additional error +Eg appears. Combining with stability yields the Wasserstein bound
(1 (7. (7.learm)
11) S Cstab (T + tE tg).
B.2 PROOF OF STABILITY
The goal is show that the composition of the propagators bound:
n II
R t lotst
S Me",
C(B)
where a Banach space of observables (c.g., equipped with the norm 1 + (1D- M 2 1, and w E R. This ensures that the operator splitting scheme remains stable over multiple time steps.
To ensure the stability the Lie-Trotter splitting scheme, we need to bound the operator norm of the composition of the propagators S Sic-1+ S((A) Stti+t over steps, where S(4) and (BB) are the solution operators (propagators) associated with the flows (A4) and ((B) acting on a Banach space of observables B. We define B as the space of test functions & U E equipped with the norm: || = || + ||D where || supneu ((((t) and D supneu [.D.(())lu-
((A)
Proof. We proceed by analyzing the action of the propagators ((A) Sthee and Sthe- over a single time step [t,t + h] with h r. The propagator S corresponds to the flow o (1) evolves the field f according to:
(+++D which
ds f(s) (y
g(s) = g(t), SE (t, h..
and similarly for (BB), For function to B, the propagator is defined as:
s (A) ( (())
We need bound the operator norm:
| (4)
(18(4)
(18 + (ID(S((A) (1)))++))))))))))
sup
sup
|
16

---

Published as a conference paper at ICLR 2026
B.2.1 BOUND ON THE SUP-NORM
Since ((A) (the- is flow map on 24 we have:
11tth sup (p( ( (+hee (u)) z sup (pp(')| || well
(4)
w'Eut
because the flow map simply reparametrizes the domain of yp.
B.2.2 BOUND ON THE GRADIENT NORM
The Fréchet derivative of ((A) i- is:
= Dip( D
(A) (u)) (
where D' (A)
(u) : + u is the Jacobian of the flow map. We need to bound
| (4)
(w) The flow ((A) satisfies the ODE:
Φ 1 (s) 1 (22.())) (u)),
(1) (1)
ds
The Jacobian (A) (u) (u) evolves according to the variational equation:
D4 (4) ds D (s, (2)( (2l2.(u,
d
(2)2((u)= =
By assumption (A1), (t,) = (t, 0) satisfies the Lipschitz condition:
[|o(A) u) - (A) u)llu (t, f.9) ft s Lallu ullu,
for ū in a ball of radius R. Thus, the operator norm of the Fréchet derivative Dol (1, u) is
bounded by LR. Applying Grönwall's inequality to variational equation:
D (A) thh (1)llee) S w Lads hnch
Therefore:
| sup D well
S DiPlllooe
Combining both bounds:
[|S(A] | I 1|o| + Dypllo S
+ [1D()la>) = chark [
80:
|S(A) p (E(B) S etin
Similarly, for (B) Sithet obtain:
| lC(B) S bnun
since (B) satisfies the same Lipschitz bound by (A1).
B.2.3 COMPOSITION OVER ONE STEP
For the composition (BB) St o S(A) Stpheri ||
(A) (1((B) llstheelll(By (|4theello(BB & Sehah Ebah 2Lah
17

---

Published as a conference paper at ICLR 2026
B.2.4 COMPOSITION OVER STEPS
Over 1/T steps covering [0, 1]. with T-= 1/n, the total composition is:
R-1
R-1
πI k=0
S(A) 1
aLat 22 = 22Lx11 = 2LRR
((B)
k==
Thus, we set M 1 and w nL satisfying the required bound:
R-1
π
s + o s(d)
S S
L(B)
B.3 PROOF OF CONSISTENCY
Consistency the local (difference between the exact flow and the split flow one subinterval) vanishes as h + O. For every t E [0, 1) and T,
S 24theet v)) L(1)g
(ith+
lim h-0
h
uniformly on bounded subsets of U.
Proof. Fix a bounded set K C l and u E K. We expand ((A))p(u)) (u)( using Taylor's theorem for Fréchet differentiable functions. For small h, the A-flow satisfies:
(AA) (+h+c( = hx (4) (1, u) Ra(h, u),
where the remainder (h, u) satisfies (0,000) chitzness of v).
h 0 (uniformly on K, by local Lips-
Expanding ip around u:
((((4) = :p()) + (t,) + Rahh, u)) + +(ll/u(()) - RAll).
Dividing by h:
( 1 ((())((() p(u) p(u)) (Dup(u), (t,) u)) 1 (Duplu), Ralp,u),,Rahh, u)) + [1hoolA] + RAll
(
h
As > 0, the second and third terms vanish (uniformly on K), so:
1 h
(s(()p(()) v(())) (Dup(uu) (1)) u)) uniformly on K.
Now apply the B-llow to S(A) Expanding (1B) ((A) p(u) (0(()) (u))) similarly:
(B) (th ( (1()) (u) o (() hw B3) (()) - Rs (),
with Ra 0 uniformly K. Expanding & around ((A)
(( (T(A)) (u))) (()( h (T x (B) (())()))) o(h)-
Substituting the A-step expansion:
S(A) p(u) = =(u) + h(Dw(((), e (A) ((a) u)) + t v (B) (1, w)) - o(hi)
since Dusp( (A)) (u)) Dusp(u) + o(1) (continuity of Dusp) and v (B)) (()) w(B) (1,
o(1) (contimuity of v).
Combining terms and dividing by h:
1 (s(B)g(A)(())) k(u) (Dup(u), u(t, + o(1) -c(())(u)+++ + oo))
h
Taking h 0 gives the result.
18

---

Published as a conference paper at ICLR 2026
B.4 PROOF OF CONVERGENCE
Below, we provide a complete and rigorous proof of the convergence result stated in the paragraph. The proof is constructed based on standard results from operator splitting theory for ODEs Hilbert spaces, eιror estimates for Lie- Trottert methods, and bounds on Wasserstein distances for pushfor- ward measures. It assumes the conditions from Assumptions (A1)-(A3) as defined in the document. The proof is divided into steps for clarity: first, bounding the splitting error with exact vector fields, then incorporating the learning error, and finally combining them via the triangle inequality.
B.4.1 PRELIMINARIES
Recall that the coupled evolution is governed by the flow map 4. : 2 - U. satisfying the ODE
d
di (e(w)) =o(t, (r()),
Po(u) u, u (f.g) e,,
where o(t, -) = (Vy(t, 9.9). ftt. the velocity field from the weak continuity equation Eq. 2 The true measure at time is la + #Po (pushforward of μn under 1
The Lie-Trotter splitting approximates this flow with component flows {(A)) and $ (B) where: For $ (A), d D (u) (t, ( with (t,) u) (V,(t, 0), For क [B), Φ ()( (B) (u)) with (B) u) = (0, (t. f,9)..
(u)
((B) Φ
The split flow over one time step h is + (splat) $ o (A), Over the interval [0, 1] with N 1/T
(number of steps) and step size h т, the approximate measure from exact splitting is μ N (o(owlt) #1o (pushforward under N-fold composition of a (split)).
For learned velocity fields iy and ig. the learned component flows are ((A) (governed by ((A) (0y, 0)) and (b,,) (governed by (( (0, ky)). The learned split flow is & (splst)
(1²) (B) (A)
and the corresponding learned approximate measure is μ4 (T.leam) = (á.
N #
We assume HD E P2(()) (probability measures on t with finite second moment) and use the Wasserstein-1 distance:
(1( = sup | pdu - / ōdv
11
where [ppl supu,aeu uff
(d()) [()) denotes the Lipschitz constant of pb | ulke
B.4.2 STEP 1: LOCAL SPLITTING ERROR WITH EXACT FIELDS
Under Assumption (AI) (local Lipschitz continuity and polynomial growth of v), the flows 4.,
(14), and (B) exist and are unique for all / E [0. 1] (Picard-Lindellöf theorem for ODEs in Banach spaces).
To analyze the local error (error over one step h). we use Taylor expansions of the flows around w:
For the true flow: T (u) a + holl, u) + n² 2 di v(t, u) O - For the A-component flow:
o (() u + heylA) u) + h² d
di dl(A)((, u) + 0(h3), For the B-component flow: (B) (())
u hu (t, (t, ū) + 2 dt
422 d (B) (t, ū) + 0(h3) (substitute ü = 4 (u)).
Substituting /b/. (a) into ( the split flow expands to:
(split) (u)= -u h ( + (
2 d (A)
( + (o(4)) (00)))))) + O(h3),
+
2
di
di
where (2) BB) = (B) v((4) ()(() (,) is the Lie bracket (commutator) of v)(1) and (
(with Du denoting the Fréchet derivative of v).
19

---

Published as a conference paper at ICLR 2026
Since 2(4) + (BB) and at de ((4) the split flow is:
d u( (component wise additivity), the local error of de
122 (o(A) 0 [t, u) + 0(*))
(u) (pppltt (u)
2
Under Assumption (A1), the Lie bracket (b(4), v b(B) is bounded by constant Kn > 0 on any ball BR = {u U ||u]| & R} (local boundedness follows from local Lipschitz continuity v).
Thus, the local error in the I| [lu-norm satisfies:
(1 - ((pplit) O0(")==((h))
For observable propagators (acting test functions 6)) define S/ s ((u) =((), (3) (u)). The local splitting error for propagators is:
pp(u) (
(p(0 and
(p: u) = h (s
( (pt)) p(u)) [(t)p(u) = O(h).
where [(I) - ((, ) D is the generator of the true flow. This holds for cylindrical test functions sp T with bounded Fréchet derivatives (Assumption (A2) on well-posedness of the continuity equation).
B.4.3 STEP GLOBAL SPLITTING ERROR WITH EXACT FIELDS
To extend the local error to the global error over [0, 1], consider the N-fold composition of the split flow:
(plopiit)
(eplit)
T (eplit)
(split) (N times), with N7 1.
By the Lady Windermere's ' Fan Lemma (a standard tool for operator splitting (Hairr et al. 2006)) or telescoping sum argument for splitting methods, the global error of the flow map satisfies:
Φ (u)
(ab(ontit)
N (n) S Ce"I,
where: -C depends on the commutator bound KR (from Step 1) and the time horizon L >> 0 is the local Lipschitz constant of (from Assumption (A1)).
The bound assumes trajectories stay in a bounded ball By (justified by the polynomial growth of and finite second moment of Po. which prevents trajectories from escaping to infinity).
To translate this flow error a measure error (Wasserstein-1 distance), use the definition of W1 and Lipschitz continuity of test functions:
Wi Wi (7) = W1 #o. (ap(oplit) N #1o) -
/ ((()())) - (((((ppt)) N (w) dpo(u). sup || l
By Lipschitz continuity of de ( - ((a(ooute)) N (u) 15 1 (o((ppitt) N (w)lke-
Taking the supremum over and using the flow error bound:
(1 < [|]](u) - (aD(eptu)) N
(w)lleduo(u) < C.tab
where Cstab = Ce incorporates stability constants (from Section 3.3 and the time horizon. This is Chermof--type convergence result in the weak sense, extended to W1 via Lipschitz bounds.
B.4.4 STEP 3: LEARNING ERROR IN FLOWS
Under Assumption (A3) (accuracy of learned velocity fields), the learned fields éy, in satisfy uni- form bounds on bounded sets BRi
[V, u) (t, w)|F= S San, (tv u) G u))l S Eg.
for all E [0, 1] and u E Br. Define the total learning error as Eno = E + Eg. So:
||o(t, u) o(t, u)lu S Eua VI [0, 1). a E Bu-
20

---

Published as a conference paper at ICLR 2026
Let De denote the learned flow (governed by 0): de (b()) = ö(t, (.(w)) with @o(u) the flow error (u) = P.(u) - Lr(u); its time derivative satisfies:
Define
dt ec(()) =(t, (u)) - att, (u)) utt, (t()) vtt, (((w)) + ptt, (u()) - ö(t, Pe(u))
T2
Bounding each term: Ti: By local Lipschitz continuity of (Assumption (AI)), [Tilla S
[R]le< T2: By Assumption (A3), [[allu S Stot-
Combining these, the error ODE becomes:
dt (ex(l))lxx S Enlles(u)]le + Eux.
Applying Gronwall's Inequality (with initial condition [|eo())]lu = 0): [eex())llr S Ext
For 1 (final time), this simplifies to:
[[ & (u)lle S E
Translating to a measure error (via W1), similar to Step 2:
Wi # & #o) S Eiee
B.4.5 STEP 4: TOTAL ERROR FOR LEARNED SPLITTING
We now analyze the total error of the learned splitting scheme, where the learned approximate measure is p1 (r.leam) = {5)
(onlit))
N # Here, ((splät) = + ( denotes the learned split
flow: ((A) evolves under the learned field ((4) = (Ø,, 0). and ( evolves under (( (0, 0g).
First, we bound the splitting eπrоr of the learned flow (analogous to Step 2 for exact fields). Under Assumption (A3) (learned fields are locally Lipschitz with constant Lmontel Re LR. matching the Lipschitz regularity of the true fields), the key properties of the split flow are preserved:
The Lie bracket of the learned components (ū(A), G(B)] is bounded on bounded sets (by K'R R s KR. since Lipschitz continuity implies bounded Fréchet derivatives).
The global splitting error for the learned flow ($b(tt)t) N
(relative to the learned full flow
41) is O(T), with stability constant Õstab Cstabb (the smallness of Eno ensures Lmodal and
KR R remain close to Lg and KR).
Next, we bound the error between the exact split measure fi and the learned split measure (r.leam) This error arises from the discrepancy between the exact split flow ((s(eplit) and the learned split
flow (i(sppitt Using an argument analogous Step 3 (Gronwall's inequality for split flows):
For each step & E 1.. - 1}, the error between Φ (splie) and (split) is bounded by ((EwIT) (local error), Compounding this over N = 1/T steps gives a global flow error of ( with constant Ć' (depending on Lmoael and T = 1).
Translating this to the Wasserstein-1 distance (via Lipschitz continuity of test functions and finite moments of po):
Wi (ui 14 (r.)am))) ČE
Finally, we combine the two eπrr components using the triangle inequality for the Wasserstein-1 distance:
Wi [1] (r,laam)
111) s Wi ( (r.leam)
(11
.1(") + Wi (ni),,)) .m..
21

---

Published as a conference paper at ICLR 2026
Substituting the bounds from Step 2 (Wi (u1 (7) (un) S CstahT) and the above learned split error:
W1
p1) < Cε Caaab..
Defining Cio max{Caub. C} (a constant depending the Lipschitz constants LR. Lmodul. com- mutator bounds KR, K stability constants M, w from Assumption (A3), and time horizon = 1), this simplifies to:
(e.m)) W1 [11
p11)
S Coo (T E1oi)
This completes the convergence proof. A key implicit assumpion-ttaajcctoiis remain in a bounded set BR where the uniform Lipschitz and learning error bounds apply-is justified by: The polyno- mial growth condition on v (Assumption (A1)), which prevents trajectories from escaping to infinity in finite time, The finite second moment of fio (hence finite first moment). which ensures the bulk of trajectories stay within BR for some R > 0.
For infiniteddimensioonal spaces U (e.g., function spaces), the weak formulation of the continuity equation (Assumption (A2)) ensures that all integrals and pushforward measures are well-defined (measurable and integrable), as cylindrical test functions avoid issues with infiniteddimensional norms.
C EXPLANATIONS IN THE FFM FRAMEWORK
C.I OVERVIEW OF FUNCTIONAL FLOW MATCHING
Functional Flow Matching (FFM) is a function-space generative model that generalizes the recently- introduced Flow Matching model to operate in infiniteddimensional spaces (Kerrigan et al. 2023) The key motivation for developing models in functional spaces stems from the fact that many real- world data types, such as time series, solutions to partial differential equations (PDEs), and audio signals, are inherently infinite- meinite-dimensiinaal and naturally represented as functions. Traditional gener- ative models that operate on finite-dimensional discretizations (c.g., grids) are often tied specific resolutions, making them ill ill-posed in the limit of zero discretization error and difficult to trans- fer across different discretizations. In contrast, functional approaches like FFM are discretization- invariant, allowing generation of functions that can be evaluated at any arbitrary points, ensuring consistency and flexibility for infiniteddimensional data.
FFM works by defining path of probability measures that interpolates between fixed Gaussian measure and the data distribution, then learning a vector field on the underlying space of functions that generates this path. Unlike density-based methods, FFM is purely theerrte avoiding issues with non-existent densities in infinite dimensions (e.g., no analogue of Lebesgue measure exists in Banach spaces). It enables simulation-fres free training via regression objective, making it efficient and suitable for function spaces.
Let X C d and consider real separable Hilbert space F of functions f : x - R equipped with the Borel r-algebra B(F). The goal is to sample from a data distribution v on F. FFM constructs a path of measures (He)) where po is fixed reference measure (e.g., Gaussian) and (41 v.
A time-dependent vector field [0, 1] x F F generates the flow & [0, x F F via the ODE:
(,), (()(g)))) Φ6(g)
The path is Me [dx] the pushforward of μoo
C.2 MATHEMATICAL JUSTIFICATION MARGINAL-TO--IINS CONSTRUCTION AND SPLITTING
We recall the marginalization identity established in functional flow matching (Theorem I of Ker- rigan et al. (2023)): if for each f the conditional path uf is generated by vector field of and is absolutely continuous with respect to the marginal pt, then the averaged vector field
0, (u) = 1 ouf(1)) (u) duf (1) du (f)
duc
22

---

Published as a conference paper at ICLR 2026
generates the marginal path 14 pá dul (f), (f), and the pair (6), satisfies the weak continuity equation.
In our split setting, apply this construction coordinate-wise. For almost every joint sample (f. from v, let of (f11)) be the conditional velocity for f with frozen, and let w ((91f)) be the conditional velocity for g with f frozen. Lifting these to the product space by zero- vero-padding, we define
(A) (f.)) = (ʃ (f (/
duf (f) du d (f), 0). du(lf)
and analogously
ot (f. := (o. f /w/olin (
duf (g) (o) dur(g)
dv"(g'))
Under absolute continuity and integrability, these define the A- and B- component fields, yielding generators
Ca(t) (v: (ot DDr), Cb(t) = (D.)
and hence the full Liouville operator L LA + LB. This provides principled justification for learning of and vàς separately and composing their lifted actions in the joint space.
C.3 EQUIVALENCE TO PAIRED-TRAINING SCHEMES VIA BIJETIVE REPARAMETRIZATIONS
Several applied works, such as the marginal-to-PDE framework of Zhang et al. (0025,, train con- ditional models with mixed parameterizations. For example, one may train ii using pairs (ft, 91) and of using (f), where fi, gi denote prescribed boundary or terminal states. Despite appearing different, these paired schemes are mathematically equivalent to the decoupled training described above, once we account for bijective reparametrizations between the complementary coordinates mapping 9r 91 along the frozen B B-flow).
Theorem C.1 (Equivalence under bijective coordinate reparametrizzaionn). Let Ry 9 G be measurable bijection such that for fachfrozen f, the reparametrization between the complementary coordinates in the two training schemes given by Ry (and analogously R2 for the f-coordinate), Assume Ry, R,y and their inverses are Lipschitz on bounded supports and preserve conditional path measures. Then training of with pairs (fr, is equivalent, up to reparametrization, to training with pairs (Jt, ge); similarly for vz. Consequently, inference using oi{ wu trained under either scheme yields the same joint pushforward after mapping representations appropriately.
Proof. Let V denote the representation map from scheme (A) to scheme (B), e-g. (g() 91- By assumption, V is bijective with Lipschitz inverse and pushes forward conditional measures correctly. For scheme (1) (pairs gc)) the training loss is
(0,) E, E grd [ller (g) uf (1:
Changing variables h = =g,(9) gives
[1(0)) = En #uf (1r (
Define of (1) = v/ ( Then this loss coincides with the scheme (2) objective (pairs (fe.91 = h)) with target of. Thus the two schemes are equivalent up to pullback by My The same argument applies to v.. At inference time, pushforwards commute with coordinate changes, so the joint pushforward is identical under either parametrization.
D ADDITIONAL. DETAILS FOR DATASET
Synthetical dataset. The "Easy distribution" is defined as two-component Gaussian mixture. One component (weight 0.6) is centered at the specified mean, while the other (weight 0.4) is slightly shifted. The covariance matrices capture different correlation structures, with one showing a positive correlation of 0.8 and the other negative correlation of -0.42, along with distinct variances across
23

---

Published as a conference paper at ICLR 2026
dimensions. The "Complex distribution" is constructed from multiple equally weighted Gaussian components arranged along a circle. These components feature ingl-ddeenndet correlations, vary- ing standard deviations, and a small isotropic Gaussian perturbation that enhances multimodality.
To construct datasets that preserve conditional information, we avoid direct joint sampling and in- stead generate two complementary datasets: one focusing p(y[æ) and the other on p(æ) lyly. ) A portion of samples is drawn directly from the mixture to preserve global structure. For the p(yle (z) dataset, x drawn uniformly from [-3, 33, and y sampled from the conditional distribution ply (r). Conversely, for the p(zly) dataset, y is drawn from scaled Beta distribution over [-2, 21, and x is sampled from p(r y). To further perturb the marginals, we introduce a mixing parameter that blends the original x marginal with a uniform distribution and the original y marginal with a scaled Beta distribution, while keeping the conditionals intact. The conditionals for Gaussian distributions are computed exactly and those for Gaussian mixtures are obtained via mixture- conture-conditionnng.
Turek-Hron. This simulation models FSI scenario using the CFD solver (Lilypad (Weymouth & Yue (0TT)) where flexible Bernoulli beam is placed downstream of a circular cylinder a two-dimensional viscous flow. The system features Reynolds number of Re=1000. Let us make D denotes the diameter of cylinder. The spatially computational domain is 20D x 14D. The flexible beam has length of 2D and thickness of 0.06/), characterized by an elastic modulus of approxi- mately 336,000 and density of 10. The simulation captures the coupled dynamics between the fluid flow and structural deformation over 10 time units episode with fixed timestep of 0.001. Key outputs include velocity fields, pressure distributions, and beam displacement data saved at each time step for experiments.
Normally, the original solver addresses a FSI problem, with one FSI trajectory computed as ground truth. To generate decoupled data, we implemented two solver modifications. First, for solid-to- fluid unidirectional coupling, we compute flow fields using only a single square cylinder (the side length is D) flow, then inversely calculate beam forces from this fluid data to control the beam, and finally compute resulting flow fields induced by its motion. Second, for fluid-to-solid unidirectional coupling, we generate flow fields using a 0.5D cylinder at each timestep to inversely compute beam positions. Both datasets contain no bidirectional coupling data, and crucially use flow conditions (square cylinder and reduced- diametrr cylinder) completely distinct from the test set to prevent information leakage.
This scenario has extensive applications across diverse fields involving flexible structures in fluid en- vironments. The methodology proves particularly valuable in aquatic ecosystem research, where it accurately models the dynamic behavior of submerged vegetation such as kelp forests and seagrass beds under varying current conditions, enabling scientists to understand their role in coastal pro- tection and marine biodiversity conservation. In biomedical engineering, the framework facilitates critical hemodynamies studies by simulating blood flow interactions with deformable arterial walls, helping predict aneurysm rupture risks and optimize cardiovascular device designs like stents and artificial heart valves. The renewable energy sector benefits from applications in wave energy har- vesting, where flexible oscillating elements mimic the motion of aquatic plants to efficiently capture ocean wave energy.
Double Cylinder. This simulation models dual-cylinder fluid-structure interaction system where a fixed upstream cylinder is positioned ahead of a freely oscillating downstream cylinder that can vibrate in the y-direction. The flow is characterized by a Reynolds number Re=1000. Let us also make the diameter of cylinder D. The computational domain spans 8D x 8D lengths (128 x 128 grid points). Both cylinders have the same diameter D. The fixed upstream cylinder is positioned at (1.5D, 4D) while the oscillating downstream cylinder is located at (4D, 4D), creating center-to- center spacing of 2.5/. The downstream cylinder exhibits a mass ratio 15.0 (relative to displaced fluid mass), dimensionless natural frequency is 0.01, and damping ratio is 0.8. The simulation employs a time step 0.1 over a total duration of 1000 time units, capturing characteristic vortex- induced vibration phenomena. Our method of generating decoupled data is similar to the Turek- Hron.
This scenario finds critical applications in underwater structural stability and molti-body hydrody- namic systems across various engincering domains. In offshore engineering, the model accurately predicts the complex wake interference effects between multiple platform columns subjected to
24

---

Published as a conference paper at ICLR 2026
ocean currents, where the upstream cylinder's vortex shedding can trigger vortex- vortex-induced vibra- tions in the downstream cylinder, leading to structural fatigue. Additionally, the simulation sup- ports the design of underwater cable systems with multi-point anchoring configurations, deep-sea aquaculture net cage arrays where inter-structure spacing affects fish habitat stability, tidal energy converter farms requiring optimal turbine spacing to balance energy harvesting efficiency with struc- tural safety, and autonomous underwater vehicle swarms where wake interference affects formation stability and energy consumption.
It is worth noting that in our FSI experiments, when generating the structurecoondttion-nn-- decoupled data, the prescribed fluid fields had strong periodicity, resulting in negligible structural deformation. This led to a limitation: the model never observed meaningful structural motion during training, and thus could not reasonably be expected to predict structural dynamies. To address this, in all experiments we replaced the pruucure-conditionnte fluid decoupled data with fully coupled data, while still using the huiddcondition-on- -nnssuuutnn---tttuuuuueeee decoupled data (constructed as described above). Importantly, this "half-decoupled" setting is fair compromise that does not undermine the core functionality of our method nor diminish the demonstration of its key advantages.
Neulear Thermal Coupling. The NT Coupling dataset is generated using the open-source MOOSE (Multiphysies Object-Oriented Simulation Environment) framework (Icenhour et al.) 2018 The governing equations for the coupled evolution of the neutron flux (Ø), fuel temper- ature (T.), and fluid temperature, velocity, and pressure (Ty, L. p) are provided in the following equations. The physical fields are discretized on separate meshes: 64 x 8 for the fuel temperature, 64 x 12 for the fluid variables, and 64 x 20 for the neutron field. Neutron physics and fuel temper- ature are solved using the Finite Element Method (FEM), while the fluid domain is computed using the Finite Volume Method (FVM). All fields subsequently interpolated onto common grid to ensure spatial alignment. An adaptive time-stepping scheme is employed, and 16 output frames are recorded to capture the transient evolution of the system.
1 OΦ(2, v) DAØ + (bE, - Ex x 10, L. - Ls}, y E 10, Lyl. E [0,5].
0t
(6)
v(0,v.t) = t(.
C(L, + Lg. = 61,, 0) = 6(1x, Ly t) = O.
OT (T,(H,
pCp
V (k. VT,) Ao..
x xE 10, Lal, y E 10, Lyl, E [0,5].
Øt
0T,
aT,
(x. 0,t) Qy
(r. Lyg t) 0,
Oy
T3 T, V, = TT y.t)
v 0,
x E (b., L, + T.p.. y E 10, Lyl, L [0.5].
+ T. Vū = -p 2 u f. E L. + Lg}. y E (0, Ly). E [0.5].
0t
(8)
OT,
+ t +z.VI,) = k, V2T, x [L. L. + Lyl. E 10, Lyl. te [0,5].
PCp
Ot
oT,
ar.
kj dx (Le, y) = ke Ox LL, y,
We generate the data following the pre-iteration method (Zhang et al. 2025 ensuring the resulting dataset preserves the intrinsic multiphysics interactions. The sequence begins by holding the fluid variables and fuel temperature constant while solving for the neutron field. The resulting neutron distribution, together with the assumed constant fluid temperature, is then used update the fuel
25

---

Published as a conference paper at ICLR 2026
Fuel
Fluid
v
Neutron
Figure 4: The geometric configuration and coupling mechanism of NT Coupling setting.
temperature field. The process continues by computing the fluid fields using the partially coupled states, followed by final fully coupled pass over the neutron and fuel temperature equations. For simplicity the coupling process, we use spatially distributed fuel temperature fields as inputs in place of the explicit heat flux.
E ADDITIONAL VISUALIZATION RESULTS
In this section, present additional visualizations that omitted from the main paper due space limitations. Figure 5 shows the coupled simulation results for the same sample as in Fig- ure 3 but using CNO as the backbone. Figure provides analogous visualizations using FNO* as the backbone, on sample from the Double-Cylinde dataset. In both cases regarddesss of back- bone or dataset, our proposed GenCP consistently delivers visually superior results compared the baselines.
Furthermore, illustrate the difference between decoupled data used for training and coupled data used as inferring, Figure/Jreports predictions of four physical fields under three settings: (i) training and inference both on decoupled data, (ii) training on decoupled data but inferring on coupled data using GenCP, and (iii) joint training and inference on coupled data. The results reveal that the coupled fields differ significantly in pattern from the decoupled ones, especially near structures. Nevertheless, our GenCP successfully learns the underlying physical interactions from decoupled data and incorporates them into the mot-hha generative process, thereby producing predictions that closely approximate the truly coupled solutions.
Figure presents the visualized results of coupled inference with decoupled learning, using different paradigms NT Coupling setting.
26

---

Published as a conference paper at ICLR 2026
Surrogate-CN CNO
M2PDE-CNO
GonCP-CNNO
Figure 5: The visualizations of results from different paradigms using CNO backbone Turek- Hron setting.
IMPLEMENTATION DETAILSS
F.1 PSEUDOCODE FOR TRAINING
Algorithm Training decoupled conditional velocity fields.
Input: decoupled datasets D,, Dy
2: Initialize 0,, 0
F do 4: for batch (f1, 9n) ~ DD do
3: for epoch
5: Sample 1 (0, 11. of ~ nF, kg ~ Tg
6: Compute ft, St. ft and update 0, to minimize Il. 4;
7:
end for
8: for batch (fo, 91) ~ D. do
9: Sample (1, 1), 2'y ~ nF; 2/y ~ ng
10: Compute ft. Dt. år and update 0, to minimize l|.. - (, ge. t: f2)1
11:
end for
12:
end for
27

---

Published as a conference paper at ICLR 2026
Surrognte-ENOO
MAPDE-ENO
GenCP- CFN-
Figure 6: The visualizations of results from different paradigms using FNO* as backbone on Double Cylinder setting.
F.2 BACKBONE MODELS
This section provides detailed descriptions of the two backbone models used in GenCP frame- work: Fourier Neural Operator with modifications (FNO*) and Convolutional Neural Operator (CNO). Both models operate in functional spaces, making them particularly suitable for interaction problems where the solution fields are naturally represented as continuous functions over spatial domains.
# The FNO* model operates in the functional space by leveraging the Fourier transform to represent solution fields as combinations of Fourier modes, providing resolution independence, global receptive fields, and alignment with the spectral bias of neural networks toward smooth lutions typical in physical problems. Our implementation combines the Scalable Interpolant Trans- former (SiT) framework with Fourier Neural Operator components, where the core computational unit is the spectral convolution layer operating in the Fourier domain as F-1(R. (x), with R being learnable linear transformation and the model truncating high-frequency modes by only learning weights for the lowest M1 x M2 x M3 modes. The architecture incorporates compact ID FNO-based final layer that operates along the patch dimension through the sequence OutputProj(GELI /aapppPPooo/LLL///aaaaooo////// (SroomtSppettaaonon ld(z) + Conv (onv)d())))), along with adaptive layer normalization for timestep conditioning AdaLN(z, c) = (c)) LayerNorm (x) + (((c,, and both spa- tial and temporal sinusoidal positional embeddings for rectangular domains and temporal sequences.
Training the #NO* model employs AdamW optimizer with learning rate x 10-4, co- sine annealing with warm restarts, batch sizes of 8, gradient clipping with maximum norm
28

---

Published as a conference paper at ICLR 2026
Field
Field
Field
Field SDF
Decoupled
Coupled
Joint
Figure 7: Comparison of results across different training strategies: Decoupled, Coupled, and Joint.
Each column corresponds to physical field (u, v. p, and SDF).
of 1.0, Xavier uniform initialization for linear layers, and zero initialization for output lay- ers to ensure stable convergence. The model processes inputs through spatial patch embed- ding combined with positional encoding. timestep embedding via sinusoidal encoding followed by MLP transformation, and alternating spatial-temporal transformer blocks where spatial tention operates on (b x f. num -utchhs hidden_ding dim tensors while temporal attention processes (b num patches, frames, hidden dim) configurations, ultimately producing outputs through the compact FNO final layer that maintains spectral properties while reducing computational complex- ity through spectral convolutions along the patch dimension.
CNO. The CNO model operates functional space through a hicrarchical U-Net-like architec- ture that preserves the infiite-ddimensional nature of solution operators by employing anti-aliased operations, multi-scale processing, and explicit functional lifting and projection operations that map between finite-dimensional representations and infiniteddimensional function spaces. The chitecture consists of CNO blocks with filtered Leaky ReLU activations to prevent aliasing dur- ing upsamplingddowssampiing operations, residual blocks incorporating Feature-wise Linear Mod- ulation (FiLM) for timestep conditioning as FiLM(rx, t) I (1 + r(t)) + B(t), and lift- project structure where the lifting operation Lift Ru Rdua maps input functions to higher- dimensional latent space while projection Project Rdam - Bh maps back to output space. The network processes inputs through five distinct phases: lifting so = Lift(( encoding with downsampling CNO blocks and residual skip connections I;+1 = CNOBlock (P) (x;) and ResBlock-NNs (xx)) bottleneck processing Tneck ResBlock Nesa (Jadeepe)) (TTaeepe) decoding with up- sampling and optional invariant blocks r; CNOBlock (U) (Concat(I,, Sx-n):))(; SN):)). and final projection Moutput Project Concat(ro, so)).
Training the CNO model utilizes Adam optimizer with learning rate x 10-4, exponential decay scheduler with 7 0.99, batch size for 3D volumes, batch normalization without dropout for regularization, and MSE loss with gradient penalty for stability. The model integrates timestep infor- mation through sinusoidal encoding femb MLP((sin (w)t(s))( cos(wtt), ... sin(wa,;2)) cos cos(wa/2t)) and employs carefully designed anti-aliasing filters with cutoff frequency w. s/2.0001 slightly
29

---

Published as a conference paper at ICLR 2026
Neutron
Fuel
Fluid
Prp
TX
Grane the
10
1.43
UNI WOLSTTTO
Campan
Figure 8: The visualization of the final frame of predicted results from different paradigms NT Coupling setting. (Here we only represent Ty as fluid field visualization.)
below the Nyquist frequency, half-width A = 0.8s 8/2.0001. 6-tap Kaiser window filter design, and upsampling factor of 2 for LReLU operations to maintain the continuous function interpretation throughout the hierarchical processing while capturing both fine-grained local features and global structure necessary for accurate multiphysics interaction modeling.
F.3 HYPERPARAMETERS OF BACKBONE MODELS
To ensure a fair comparison and minimize bias from hyperparameter tuning, we standardized the backbone hyperparameters for each task. Minor adjustments made only accommodate the inherent modeling complexity of different methods, as detailed Tabel
Table 5: Hyperparameters of CenCP-FNO* for different datasets.
Turek Hron
Double Cylinder
NT Coupling
Hyperparameter name
fluid
structure
fluid
structure
fluid
fuel
neutron
Number of Blocks
4
4
4
FNO mode
4
4
4
4
4
4
Hidden size
256
128
256
128
256
128
128
Patch size
[2, 2]
[2, 2]
[2,
2]
[2,
2]
[2, 2]
[2, 2]
[2, 2]
Number of blocks
MLP ratios
4
4
4
4
30

---

Published as a conference paper at ICLR 2026
Table Hyperparameters of M2PDE-FNO" for different datasets.
Turek Hron
Double Cylinder
NT Coupling
Hyperparameter name
fluid
structure
fluid
structure
fluid
fuel
neutron
Number of Blocks
FNO mode
4
4
4
4
4
4
4
Hidden size
256
128
256
128
256
128
128
Patch size
[2,
[2,
[2,
[2,
[2,
2]
[2,
Number of blocks
4
4
MLP ratios
Table Hyperparameters of SurrogateFFNO* for different datasets.
Turek Hron
Double Cylinder
NT Coupling
Hyperparameter name
fluid
structure
fluid
structure
fluid
fuel
neutron
Number of Blocks
FNO mode
4
4
4
4
4
4
4
Hidden size
128
64
128
64
128
64
64
Patch size
[2,
[2,
[2,
[2,
[2,
2]
[2,
2]
[2, 2]
Number of blocks
MLP ratios
Table 8: Hyperparameters GenCP-CNO for different datasets.
Turek Hron
Double Cylinder
NT Coupling
Hyperparameter name
fluid structure
fluid
structure
fluid
fuel
neutron
Number of layers
2
2
2
Number of residual blocks per level
1
1
1
1
1
Number of residual blocks the neck
6
6
6
Channel Multiplier
16
16
16
16
16
16
16
Size of Conv Kernel
3
3
3
3
3
latent dimension in the lifting layer
64
64
64
64
64
64
64
Table Hyperparameters of M2PDE-CNO for different datasets.
Turek Hron
Double Cylinder
NT Coupling
Hyperparameter name
fluid
structure
fluid
structure
fluid
fuel
neutron
Number of layers
2
2
2
Number of residual blocks level
1
1
I
1
1
Number of residual blocks the neck
6
6
6
Channel Multiplier
16
16
16
16
16
16
16
Size of Conv Kernel
3
3
3
3
3
latent dimension in the lifting layer
64
64
64
64
64
64
64
31

---

Published as a conference paper at ICLR 2026
Table 10: Hyperparameters of SurrogateCNNO for different datasets.
Turek Hron
Double Cylinder
NT Coupling
Hyperparameter name
fluid
structure
fluid
structure
fluid
fuel
neutron
Number of layers
2
2
2
Number of residual blocks per level
1
1
1
1
1
Number of residual blocks the neck
6
6
6
Channel Multiplier
8
8
8
8
8
8
8
Size of Conv Kernel
3
3
3
3
3
latent dimension in the lifting layer
64
64
64
64
64
64
64
F.4 INFERENCE SCHEME OF BASELINE METHODS
Surrogate- based scheme The surrogate model in this article refers to neural network that di- rectly maps system inputs to outputs one step, unlike diffusion models that generate outputs step-by-step through denoising. For multiphysics problems, surrogate models are trained separately for each physical field using decoupled data and obtain coupled solutions through iteration.
Our work follows the convention from M2PDE (Zhang al. 2025) for implementing the surrogate model combination algorithm
Picard-iteration approach. The surrogate models' combination algorithm are identical, as demonstrated in Alg 3. The relaxation factor is set to 0.5 by default, and the tolerance émax is set to by default, effectively running the iteration for the maximum number of steps without early stopping based convergence.
Algorithm 3 Surrogate Model Combination Algorithm for Multiphysics Simulation.
Require: Compositional set of surrogate model (241, i). 2.2.... outer inputs C, maxi- number of iterations M, tolerance maxx relaxation factor a.
1: Initialize constant fields 2;, = 1,...,
0
2: while < .M and Emax do
3: m = m +
4: for i 1, 1...., . NN do
5:
&
6:
2; ac 2 C*) + 2)
7:
end for
8:
<= [1; - 2,); i 1 N
9:
end while
Ensure: 21, :..,.., NN
The notation used in Algorithm 3 differs from this paper's notations and is briefly explained as fol- lows: e represents the trained surrogate model, such neural network like FNO, for the i-th physical field; Zys denotes the predictions from all other physical fields used as conditions; C' rep- resents the outer system inputs, such as initial conditions boundary values; 2; are the predicted fields for each physical process; 2; is temporary copy of the previous prediction; a is the relaxation factor, scalar between 0 and M is the maximum number of iterations; and 6max is the conver- gence tolerance threshold. The L1 norm (L1) computes the absolute difference convergence metric, which can be adapted to other norms like L2 needed.
The algorithm takes trained single-field surrogate models, one per physical process and trained on decoupled data, and iteratively combines them to approximate coupled multiphysics solution. starts with initial constant fields 2,, such uniform value like 0.5 or problem-specife priors, along with outer inputs C. In each itcration, it updates each field's prediction conditioned on the others and applies moving average for stabilization.. The loop ends when the maximum iterations M are reached or the convergence metric c, based on the LI norm prediction changes, falls below 6max-
The moving average with relaxation factor prevents divergence instability in iterative updates, particularly in strongly coupled systems where direct updates (a = 1) might oscillate or fail to converge. This approach resembles under-relaxation in traditional numerical methods like Picard
32

---

Published as a conference paper at ICLR 2026
iteration, blending new predictions with previous ones weighted by and to ensure smoother convergence.
M2PDE scheme The M2PDE framework composes multiple diffusion models, each trained on a single physical process using decoupled data, to simulate coupled multiphysics systems through iterative conditional sampling during inference. Unlike surrogate models that directly map inputs to outputs, M2PDE generates solutions via reverse diffusion process, implicitly handling interactions between fields.
Algorithm 4 Inference of Multiphysics Simulation with the Composition of Diffusion Models.
Require: Compositional set of diffusion models c244 C,t), i ..,,...,NN outer inputs C, number of reverse diffusion steps T, number of outer iterations K.
for k L...,, do ~ 10,
2:
3: for = T ....: do
4: for i 1 N do
5:
: - ( C,t)
Estimate clean sample for component
6: end for
Compose 2 = ($0... ( 8: 2' Po (2'100,1)
7:
Perform conditional sampling step
9: end for
10:
2(k) 0
11: end for
Ensure: 21 (K)
K)
The composition algorithm is identical, as demonstrated in Alg. 4. The number of outer iterations K set to 2 by default. As introduced in the surrogate model combination algorithm, M2PDE schema shares the same notations: e) represents the trained diffusion model (score network) for the i-th physical field; 244 denotes the current noisy estimates from all other physical fields at diffusion timestep f; C' represents the outer system inputs, such as initial conditions or boundary values; 2' the noisy state at timestep t; # is the estimated clean (denoised) sample for field i: T is the total number of reverse diffusion steps; K is the number of outer iterations; and denotes the conditional sampling function in the diffusion process.
The algorithm takes trained single-field diffusion models, one per physical process and trained on decoupled data, and iteratively composes them during reverse diffusion to generate coupled mul- tiphysics solution. It starts by sampling initial noise 2r from standard normal distribution, along with outer inputs C. In each outer iteration k, it performs a full reverse diffusion process: for cach timestep t from 7' to I, it estimates denoised samples for each field conditioned on the others' cur- rent states, composes them, and advances the sampling step. The outer loop refines the solution over K iterations, with the final output being the denoised fields after the last iteration.
The outer iterations (K) enable progressive refinement of field interactions, improving accuracy in strongly coupled systems. Conditional sampling ensures that each field's generation accounts for others, implicitly capturing multiphysics couplings without explicit iteration during training.
G LIMITATIONS AND FUTURE WORK
GenCP has been evaluated on two-dimensional synthetic distribution and multiphysical problems, without yet being applied to more complex systems like 3D or real-world engineering problems. Additionally, has not considered irregular geometric boundaries, which is another notable limi- tation. These limitations arise because we primarily focus on introducing our core paradigm, and are further constrained by limited data availability, which also opens up exciting and broad avenues for future research. Besides, we have not yet explored incorporating physical priors into our frame- work. Such integration may further improve our performance on out-of- distribution data, thereby enhancing both the accuracy and robustness of coupled predictions after training on decoupled data. Looking ahead, we also plan to leverage GenCP to unlock the value of abundant historical single- physics data and integrate ioss-institutional datasets for learning complex coupled physics.
33


Reviews:
[{"review_id": "YPCvcUYdnZ", "reviewer": "Reviewer_k4vA", "summary": "The submission proposed a flow-based model to handle multiphysics simulation. The idea is simply: training two flow models independently but using another variables as the condition. It reduces the training complexity. I appreciate the method, but not satisfied with current experiments part.", "questions": "Line 202, Line 207, I believe in the flow model, the vf = f1-f0, not f1-zf, the same for vg.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "I think the proposed method is easy to follow and the result is good. Also, I think the problem the authors are trying to solve is meaningful.", "weaknesses": "1. I think the experiments are not enough. It all concentrates on fluid dynamics, but it is better to combine fluid and rigid objects. It is convenient to treat the whole fluid system as a domain.\n\n2. Also for the experiments part, I feel like the submission should add a baseline: training a single v(x_t, y_t) for all variables (coupled data), to prove the decouple the data is beneficial.", "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": "hfmYpTYkrH", "reviewer": "Reviewer_sgcu", "summary": "This paper proposes a novel and elegant generative paradigm—GenCP—for coupled multiphysics simulations. By formulating the coupled physics modeling problem as a probabilistic modeling problem, the authors' key innovation lies in combining probability density evolution from generative modeling with iterative multiphysics coupling. This enables training on decoupled data and inference of coupled physics processes during sampling. The authors evaluate their paradigm on a synthetic dataset and two challenging fluid-structure interaction scenarios to demonstrate GenCP's fundamental insights and superior application performance. This is an article about applying AI to engineering physics.", "questions": "Have you considered a larger-scale experiment?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Combining probability density evolution with iterative multiphysics coupling in generative modeling is a novel approach, and the problem setting is very clear.", "weaknesses": "The theoretical part of this paper is excellent, but I think the experimental part could be scaled up further.", "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": "3LHiGRkNHU", "reviewer": "Reviewer_Uo4f", "summary": "The authors introduces GenCP, a generative modeling framework for simulating coupled multiphysics systems with applications to fluid-structure interaction. Their core contribution is training on decoupled data which does not require coupled simulation data during the training process. During inference,however, the  fully coupled system can be predicted.", "questions": "See weaknesses and:\n- For the two FSI examples, the authors just report single error values but the Reviewer would have expected a distribution due to the probabilistic approach. Are these posterior mean values ?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- According to the Reviewer, the general idea presented in the paper of integrating operator-splitting theory with generative flow models is novel. Using decoupled training and coupled inference later to reduce the cost of training data could potentially be very impactful as many AI4Science ML applications deal with the issue of spending significant compute on training data generation.\n\n- The authors present theoretical guarantees for their framework. The method is not just a heuristic but directly connected to operator splitting theory used in numerical schemes.", "weaknesses": "Despite the core idea being novel, the Reviewer is unable to recommend the paper for acceptance at ICLR due to the following major issues:\n- One of the central motivation of the paper is that the training data generation using decoupled systems only is cheaper but this was unfortunately never quantified. The reviewer recommends evaluating the cost of training data generation for coupled and decoupled generation.\n- The comparisons to baselines should also include frameworks trained on the fully coupled data to estimate the change in accuracy by using decoupled data only.\n- The authors only apply their method to, in the setting of FSI, small problems. However, as stated in the Appendix for one of their examples the authors do not use decoupled data only but coupled training data as the decoupled data contained no meaningful structural motion. According to the Reviewer, this undermines one of the central statement of the paper. If already for a small example, using decoupled training data only is insufficient, the Reviewer would expect for more complex system with more non-linear coupling the same issue. This would significantly limit the applicability of the new approach.\n- The obtained relative errors on the coupled data test set seem relatively high and would be too high for most engineering applications. The Reviewer acknowledges that the authors perform better than reported baselines but the gap between results on the decoupled validation set and the coupled test set seems to high and would discourage currently training on decoupled data in the Reviewers opinion. Moreover, the Reviewer encourages the authors to add the aforementioned baselines on a fully coupled dataset.\n- The authors use convolutional building blocks and do therefore require a regular grid. For FSI which often involves complex boundaries this is a significant limitation. In Figure 6, the cylinder is not round but the discretization artifacts are clearly visible. The authors should comment on this limitation.", "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": "aT4RuJIcwv", "reviewer": "Reviewer_CPwt", "summary": "This paper presents GenCP, a novel paradigm for multiphysics simulation that learns coupled physics from decoupled training data. \n\n**Key Innovation:** The elegant combination of flow matching's ODE formulation with operator splitting theory. The method trains separate flow matching models on decoupled data (learning $p(f|g)$ and $p(g|f)$), then composes them via Lie-Trotter splitting during inference to generate coupled solutions $p(f,g)$.\n\n**Theoretical Contribution:** Establishes rigorous convergence guarantees by connecting probability density evolution with operator splitting.\n\n**Significance:** This work bridges generative models with classical numerical methods, providing both theoretical rigor and practical efficiency for solving expensive coupled simulation problems.", "questions": "- The paper assumes decoupled data is easier to obtain (line 47), but your experimental protocol suggests otherwise. You still need to run solvers to generate it. Would it be more accurate to position GenCP's value as enabling the \"reuse\" of historical single-physics data or \"integration\" of cross-institutional datasets, rather than claiming acquisition cost advantages in general?\n\n- Experimental Design: The paper acknowledges in Appendix D that structure-conditioned-on-fluid training uses coupled rather than truly decoupled data. Could the authors provide ablation studies comparing: (1) fully decoupled training for both directions, (2) the current \"half-decoupled\" setting, and (3) fully coupled training? This would help quantify the performance trade-offs and clarify whether truly decoupled training is feasible for both directions, which is central to the paper's claims.\n\n- The success of GenCP relies on the assumption that decoupled datasets $\\mathcal{D}_f$ and $\\mathcal{D}_g$ contain sufficient information to recover the joint distribution. In your experiments, this is guaranteed by construction, since both datasets are derived from the same coupled solver with consistent physical parameters. Can you characterize the necessary and sufficient conditions on decoupled data for GenCP to work? Specifically:\n  - Must $\\mathcal{D}_f$ and $\\mathcal{D}_g$ share the same physical parameters?\n  - What level of parameter mismatch is tolerable between the two datasets?\n  - Are there failure modes where learning $p(f|g)$ and $p(g|f)$ separately is provably insufficient to reconstruct the joint distribution $p(f,g)$?\n  - How can practitioners assess the \"coupling information completeness\" of their real-world decoupled datasets before committing to training?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "- **Intuitive decomposition with strong empirical validation.** The core idea of learning coupled physics through separate conditional models is elegant and practical. The paper convincingly demonstrates this works: models trained only on decoupled data successfully capture strongly coupled dynamics, which surrogate baselines completely fail to model despite similar numerical errors. This validates that the probabilistic formulation genuinely captures coupling mechanisms.\n\n- **Rigorous theoretical foundation bridging generative models and numerical analysis.** The connection between flow matching's velocity field decomposition ($v = v^{(f)} + v^{(g)}$) and Lie-Trotter operator splitting is non-trivial and principled. The formal convergence analysis via Wasserstein-1 distance provides error controllability guarantees rare in neural PDE solvers, with explicit tracking of splitting error ($\\tau$) and learning error ($\\varepsilon_f, \\varepsilon_g$).\n\n- **Careful decoupled data construction avoiding information leakage.** The experimental design modifies boundary conditions to ensure zero coupling information in training data while maintaining physical plausibility. This proves the method learns underlying conditional physics rather than memorizing coupling patterns.\n\n- **Insightful analysis of baseline failure modes.** Beyond reporting metrics, the paper explains why deterministic surrogate methods fail: they miss oscillatory dynamics despite low numerical errors because they cannot handle mode errors and stochastic behavior. This clarifies that GenCP's probabilistic modeling addresses fundamental limitations of deterministic learning for strongly coupled systems.", "weaknesses": "- **Notation density limits accessibility.** The formulation rapidly introduces weak continuity equations and Lie-Trotter splitting without sufficient intuitive scaffolding. Key notations need clearer explanation before formal treatment.\n    \n- **Incomplete Decoupling.** The experiments acknowledge (Appendix D) that one training direction uses coupled rather than truly decoupled data due to negligible structural deformation in genuinely decoupled scenarios. This \"half-decoupled\" compromise contradicts the paper's emphasis on training exclusively from decoupled data, suggesting the method still requires partial access to expensive coupled simulations in practice.", "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)