Multi-Fidelity Flow Matching: Cascaded Refinement of PDE Solutions Sipeng Chen
Department of Computer Science, Florida State University
arXiv:2605.16118v1 [cs.LG] 15 May 2026
Junliang Liu [email protected] Department of Materials Science and Engineering, Florida State University Hewei Tang [email protected] Hildebrand Department of Petroleum and Geosystems Engineering, The University of Texas at Austin Shibo Li
Department of Computer Science, Florida State University
Abstract The source distribution in conditional flow matching is a design parameter that can be calibrated to data, not a default isotropic prior. We exploit this in Multi-Fidelity Flow Matching (MFFM), a cascade refinement framework for parametric PDE solutions: the source is calibrated to the empirical low-to-high-fidelity residual scale with local Gaussian-blur correlation, and the velocity network is conditioned on the low-fidelity solution. Conditioning makes the residual refinement problem substantially easier than unconditional field generation, while residual-calibrated source noise improves the flow-matching training geometry. A multi-resolution cascade applies the same construction independently between adjacent fidelities. After level-wise flow-matching pretraining, we fine-tune the composed cascade end-to-end with a deterministic one-step rollout, which makes one velocity evaluation per cascade level the optimized operating point at inference. The result is a learned analog of multigrid refinement that reaches the finest grid in L deterministic network evaluations per query. We validate MFFM on eight benchmarks: two super-resolution problems and six spatiotemporal forecasting tasks from PDEBench, The Well, and the FNO Navier–Stokes dataset.
1. Introduction Solving a parametric PDE from its coefficients, initial conditions, or boundary conditions remains costly even with modern neural surrogates. The operator ξ 7→ uHF is high-dimensional, geometry-sensitive, and unstable on parameter regimes far from training. Refining an alreadycomputed approximate solution is a different and strictly easier task: the residual carries far less variance than the field itself. A coarse-mesh solver, a reduced-order model, or a neural-operator surrogate already captures the macroscopic structure of the high-fidelity field; the residual δ = uHF − uLF then has smaller variance and is concentrated in the directions where the cheap solver loses resolution. For typical parametric PDEs the LF and HF solutions are strongly correlated across most of the domain. Classical multi-fidelity methods (Kennedy and O’Hagan, 2000; Perdikaris et al., 2017; Peherstorfer et al., 2018) and recent neural 1
multi-fidelity surrogates (Lu et al., 2022; Howard et al., 2023; Li et al., 2022b, 2024) have long exploited this regime. Multi-Fidelity Flow Matching (MFFM) sits downstream of operator-learning surrogates rather than competing with them. A neural operator absorbs the parametric variability of the solution manifold; MFFM consumes its output, or any low-fidelity solver’s output, on a finer grid and produces a high-fidelity refinement through a flow-matching cascade. In practice we expect MFFM to follow a learned LF surrogate: the operator handles the parameter-to-solution map at its native resolution, and the cascade refines to the high-fidelity grid. Figure 1 previews the geometric distinction between standard flow matching and the adapted-source variant. We pose the LF-to-HF refinement as conditional residual flow matching. Generative approaches to PDE-related problems (Huang et al., 2024; Lippe et al., 2023; Bastek et al., 2025; Shu et al., 2023; Baldan et al., 2026; Hou et al., 2026) typically take the source distribution as uninformative noise, so the model must recover the macroscopic component of uHF at every ODE step and pays the same transport cost as single-fidelity generation. We instead calibrate the source to empirical residual statistics: in the main implementation, iid Gaussian noise is locally correlated by Gaussian blur, normalized samplewise, and scaled by the empirical per-coordinate residual standard deviation. The velocity network is conditioned on uLF at every layer. The two design choices compose: (i) the residual-calibrated source reduces the marginal scale gap between source and target and better matches the local structure of PDE residuals; (ii) conditioning on uLF makes the residual refinement target substantially more concentrated than the full high-fidelity field. After level-wise flow-matching pretraining, we fine-tune the composed cascade end-to-end using a deterministic one-step rollout. Thus one midpoint velocity evaluation per cascade level is the optimized operating point used for all main predictions. Stacking the same construction across nested resolutions G0 ⊂ G1 ⊂ · · · ⊂ GL produces a learned analog of multigrid refinement that reaches the finest grid in L deterministic network forward passes per query, with per-level networks correspondingly smaller because each handles a narrower band of residual frequencies.
2. Background The classical multi-fidelity treatment of Kennedy and O’Hagan (2000) models the high-fidelity output as a scaled low-fidelity output plus a Gaussian-process discrepancy, and motivates the additive decomposition we adopt. Nonlinear couplings (Perdikaris et al., 2017) replace the scaling with a learned map; Peherstorfer et al. (2018) surveys two decades of multi-fidelity methodology across uncertainty propagation, optimization, and inverse problems. An earlier surrogate-modelling thread (Forrester et al., 2007) formalized cheap-to-expensive corrections in design optimization. Neural realizations of multi-fidelity correction include multi-fidelity DeepONets (Lu et al., 2022; Howard et al., 2023), infinite-fidelity coregionalization (Li et al., 2022b), and active-learning-driven multi-resolution operator architectures (Li et al., 2024, 2022a, 2020), all built on the broader neural-operator framework surveyed by Kovachki et al. (2023). These approaches model the LF-to-HF correction as a deterministic regressor; MFFM instead uses residual flow matching with a source distribution calibrated to empirical LF-to-HF residual statistics, followed by deterministic cascade fine-tuning for point prediction. 2
Figure 1: Standard flow matching (left) transports an isotropic source N (0, I) to the parameter-induced distribution uHF along curving ODE trajectories, so the model must absorb the full parametric variability of the target. MFFM (right) replaces the source with one calibrated to the LF→HF residual and conditions on the low-fidelity solution, making the residual refinement substantially easier than unconditional field generation. The full method composes L such per-level refinements across nested resolutions as a cascade (bottom; Sec. 3.4); after end-to-end fine-tuning (Sec. 3.5) one velocity evaluation per cascade level is the optimized operating point at inference.
Conditional flow matching (Lipman et al., 2023), a continuous-time generative-modelling framework descended from neural ODEs (Chen et al., 2018) and the diffusion family (Ho et al., 2020; Song et al., 2021), trains a velocity network by regressing onto closed-form per-sample conditional velocities, avoiding ODE simulation at training. Under the linear interpolation path xt = t x1 + (1 − t) x0 , the per-sample conditional velocity reduces to the constant x1 − x0 . Rectified flow (Liu et al., 2023) derives the same path from a trajectorystraightening argument and can support low-NFE inference when the source-target coupling is well learned. Minibatch optimal-transport coupling (Tong et al., 2024) further reduces regression variance. Across these variants the source is usually treated as a fixed reference distribution. Bridge methods (De Bortoli et al., 2021; Shi et al., 2023) and stochastic interpolants (Albergo et al., 2025) establish the broader formal framework within which non-noise sources are admissible.
Cascaded diffusion (Ho et al., 2022; Saharia et al., 2022) composes denoising samplers across resolutions, and multi-level Monte Carlo (Giles, 2015) achieves an analogous variance reduction in the deterministic regime; both inform the cascade design here. Image-domain cascades typically retain a multi-step sampler at every level, whereas MFFM uses adjacentfidelity residual flow matching as pretraining and then operates as a one-forward-pass-per-level deterministic refiner after end-to-end cascade fine-tuning. 3
3. Multi-fidelity flow matching 3.1 Problem formulation Let Ξ ⊂ Rp index a family of parametric PDEs Lξ [u] = fξ on a fixed domain Ω ⊂ Rd . A low-fidelity solver SLF produces uLF = SLF (ξ) at cost cL on a coarse grid, and a high-fidelity solver SHF produces uHF = SHF (ξ) at cost cH ≫ cL on a fine grid. Our primary setting is spatiotemporal PDE forecasting: uLF stacks T past frames at the coarse grid into channels, uHF stacks the corresponding T future frames at the fine grid, and refinement runs on the HF grid after upsampling uLF . This is the PDEBench setup of Section 5. The construction also covers spatial or space–time super-resolution (T = 1 for spatial fields), and we include Darcy and Burgers as such validations. (i) (i) Training assumes a paired dataset D = {(uLF , uHF )}N i=1 obtained by running both solvers or constructing paired LF–HF fields on identical PDE instances. The parameter ξ generates the data but is not provided to the model. The setting differs from operator learning in two consequential ways. Operator surrogates approximate ξ 7→ uHF from ξ alone and must absorb every source of parametric variability without a solution-side prior. MFFM operates downstream of SLF and exploits uLF as an informative summary of uHF ; the map we learn is the conditional refinement uLF 7→ uHF , which removes the dominant macroscopic component of the parametric variation. The two settings compose: a learned operator (FNO, DeepONet, multi-fidelity DeepONet, COMPOL (Sun et al., 2025; Li et al., 2024)) can play the role of SLF , and the resulting refinement provides high-fidelity predictions without re-running the expensive HF solver. We pose the refinement as conditional residual flow matching, learning a velocity model on the residual δ = uHF − uLF conditioned on uLF . 3.2 Adapted-source flow matching The construction is best understood by contrast. Standard flow matching transports samples from an uninformative source distribution to the data and pays a transport cost proportional to the geometric gap between the two. In the multi-fidelity setting we instead draw the source from empirical residual statistics, so the source and target share the same residual scale from the outset; the network learns the conditional coupling between source and target rather than rescaling between distributions of incompatible scale. We work directly on the residual δ = uHF − uLF ∈ Rn and define a flow-matching generator with a calibrated residual source. In the main implementation, the source is a locally correlated Gaussian-blur residual source: z ∼ N (0, I),
z̃ =
Kτ ∗ z d τ ∗ z) + ϵnum Std(K
,
ε = σδ ⊙ z̃,
d D [δj ], σj2 = Var
(1)
where Kτ is a depthwise Gaussian blur kernel, ∗ denotes spatial convolution, and the standardization is applied samplewise. This construction preserves the empirical per-coordinate residual scale while introducing local spatial correlation, matching the structure of PDE residual fields more closely than pixelwise independent noise. The calibrated diagonal source ε ∼ N (0, diag(σδ2 )) is included as an ablation. Following Lipman et al. (2023), we adopt the linear interpolation conditional path δt = tδ +(1 − t)ε, t ∈ [0, 1], which has constant per-sample conditional velocity v(δt , t | δ, ε) = δ − ε. 4
The MFFM regression objective JMFFM (θ) = Et∼U [0,1] E(uLF ,δ)∼D Eε∼qδ ∥vθ (δt , t, uLF ) − (δ − ε)∥2
(2)
exposes uLF to the velocity network at every layer rather than through a single conditioning bottleneck, so vθ can learn the conditional refinement structure rather than the unconditional residual marginal. Here qδ denotes the calibrated residual source in Eq. (1). At inference time, one can sample ε ∼ qδ and integrate the inference ODE dδt = vθ (δt , t, u∗LF ), dt
δ0 = ε,
(3)
from t = 0 to t = 1, returning uHF ˆ = u∗LF + δ1 . Drawing K independent source samples produces an optional stochastic ensemble. For the deterministic point-estimate results in Table 1, we use the end-to-end fine-tuned deterministic cascade described in Section 3.4; stochastic source sampling is not used for the main accuracy results. The choice of source in Eq. (1) is a normalization and correlation-design decision (an analog of the noise-schedule design space catalogued for diffusion models by Karras et al. (2022)), not the mechanism that makes one-step refinement effective. Matching the empirical residual scale aligns the marginal magnitude of δt with both endpoints throughout the path; adding Gaussian-blur correlation further matches the local smoothness of PDE residuals. Under the diagonal source, the transport scale satisfies E ∥δ − ε∥2 = 2tr(Σδ ) with Σδ = diag(σδ2 ); on Darcy at 1282 this is roughly one-tenth of the corresponding Θ(n) value an isotropic N (0, I) source would induce on unit-normalized fields (Appendix A.3). The mechanism making one-step residual refinement effective is conditioning on uLF , formalized in Section 3.3. 3.3 Single-step inference under deterministic conditioning By default the inference ODE (3) may require many integration steps. The training loss in Eq. (2) drives vθ toward the marginal Bayes-optimal velocity v ∗ (z, t) = E[δ − ε | δt = z, t], which is generally t-dependent: for two independent zero-mean Gaussians with equal covariance, v ∗ (z, t) = (2t − 1)/(2t2 − 2t + 1) z. Universal single-step claims for conditional flow matching are therefore false. What MFFM uses is the conditional regressor under uLF , which becomes simple for a structural reason specific to the multi-fidelity setting. One-step inference under deterministic conditioning. Condition on uLF and assume that the LF solver is injective on the parameter space, so that uLF uniquely determines ξ and, since SHF is deterministic, also determines uHF . The conditional law of δ given uLF is then the Dirac measure concentrated at the deterministic refinement target g(uLF ) := uHF − uLF . From this: δt = t g(uLF ) + (1 − t) ε (4) is a deterministic function of ε given the conditioning, and the source is recoverable from the input via ε = (δt − t g)/(1 − t) for t < 1. Substituting into the regression target yields the conditional Bayes-optimal regressor v ∗ (z, t | uLF ) = (g(uLF ) − z)/(1 − t), 5
(5)
which is generically t-dependent off the inference path but reduces to the constant g − ε along it. A single Euler step from δ0 = ε would then return δ1 = ε + (g − ε) = g(uLF ), the high-fidelity refinement target, if the trained network coincided with the Bayes regressor. In practice ĝθ approximates g to a residual error we measure empirically. After deterministic end-to-end fine-tuning, the one-step rollout becomes the optimized operating point rather than a coarse approximation to a many-step ODE solve; Section 5.5 shows that increasing NFE does not improve the fine-tuned predictor. We do not claim exactness; we claim the construction makes one step the appropriate operating point under (i) SHF deterministic and SLF sufficiently informative on Ξ, the standard PDE-solver setting; (ii) source independence ε ⊥ uLF , which holds by construction; and (iii) evaluation away from t = 1, avoiding the 1/(1 − t) singularity in Eq. (5). The deterministic e2e predictor uses the midpoint t = 1/2, as described in Section 3.4. Empirical confirmation. Section 5.5 evaluates the fine-tuned deterministic cascade with multiple Euler evaluations per level. Increasing NFE does not improve NRMSE, confirming that the one-step rollout is the appropriate operating point for the e2e fine-tuned predictor. 3.4 Multi-level cascade Two regimes break the single-shot argument of Section 3.3. Stochastic SHF (random forcing, sub-grid Monte Carlo, SPDEs) and non-injective SLF replace the Dirac measure with a non-degenerate conditional distribution; this regime does not arise for the deterministic numerical solvers we consider. Insufficient network capacity, by contrast, is the operative bottleneck whenever the LF-HF gap is large: a single-shot refinement asks one network to recover g across the full gap, where g has large per-coordinate variance and the conditional structure is harder to learn at fixed parameter budget. The cascade addresses this regime by splitting the bottom-to-top refinement into L adjacent-fidelity steps; each per-level gℓ has smaller variance and stronger conditioning than the bottom-to-top g, and the conditional argument above applies more cleanly per level than per single shot. Let nested grids G0 ⊂ G1 ⊂ · · · ⊂ GL form a multi-resolution hierarchy with prolongation operators Iℓ→ℓ+1 : Rnℓ → Rnℓ+1 (bilinear interpolation in our experiments). At each level the cascade applies one MFFM step: uℓ+1 = Iℓ→ℓ+1 (uℓ ) + δℓ ,
ℓ = 0, . . . , L − 1,
(6)
where δℓ is generated by an MFFM model trained on level-ℓ pairs (Iℓ→ℓ+1 (uℓ ), uℓ+1 ). Eq. (6) defines the cascade additive decomposition. Each level has its own velocity network vθℓ trained with the per-level objective Jℓ (θℓ ) = Et, εℓ , (uℓ ,δℓ ) ∥vθℓ (δℓ,t , t, ũℓ ) − (δℓ − εℓ )∥2 ,
(7)
with ũℓ = Iℓ→ℓ+1 (uℓ ), εℓ ∼ qδℓ , and δℓ,t = tδℓ + (1 − t)εℓ , where qδℓ is the per-level calibrated residual source from Eq. (1). Per-level deterministic conditioning. Training each level with the calibrated-source loss in Eq. (7), the conditioning argument of Section 3.3 transfers level by level: each network observes the prolongated state ũℓ and only needs to predict the adjacent-fidelity residual uℓ+1 − ũℓ . During level-wise pretraining, the conditioning states come from ground-truth 6
intermediate fidelities; during inference, they are produced by previous cascade levels. This creates a cascade covariate shift. We address it directly with the deterministic end-to-end fine-tuning stage below, which optimizes the composed one-step rollout used at test time. Per-level transport-cost bound. The per-level residual variance can be related to the empirical adjacent-fidelity correlation ρℓ,ℓ+1,j = Corr(uℓ,j , uℓ+1,j ). Treating adjacent fidelities 2 ≈ σ2 as having approximately equal per-coordinate variance σℓ,j ℓ+1,j , the per-coordinate residual variance satisfies 2 σδ2ℓ ,j ≈ 2(1 − ρℓ,ℓ+1,j ) σℓ+1,j , (8) mirroring the bottom-to-top calculation in Appendix A.3. Adjacent fidelities are typically more correlated than non-adjacent ones in our benchmarks, so the adjacent-level residual scale is smaller than the direct bottom-to-top residual scale, making each per-level transport task easier than the single-shot coarsest-to-finest task. Cost and parameter allocation. Total cascade inference cost is L forward passes per query. Higher-resolution levels handle progressively smaller corrections (Eq. (8) with ρℓ,ℓ+1 closer to one for finer level pairs), so the velocity network at finer levels uses fewer parameters per coordinate. The variance-allocation parallel with multi-level Monte Carlo (Giles, 2015) is direct: MLMC distributes samples adaptively across levels to minimize variance per dollar of compute, and the cascade allocates parameters across levels by the same logic. 3.5 End-to-end cascade fine-tuning The independent level-wise training objective in Eq. (7) identifies each refinement model on the training chain induced by ground-truth intermediate fidelities. At inference time, however, the conditioning state at level ℓ is produced by the previous learned refinement models, so the per-level networks are evaluated outside the support seen during their independent training. We address this cascade-induced covariate shift directly with a second training stage: after independent level-wise pretraining, we unroll the full cascade from the coarsest input u0 to the finest grid using one Euler evaluation per level and initialize the residual state at each level by zero: ûℓ+1 = Iℓ→ℓ+1 (ûℓ ) + vθℓ 0, 12 , Iℓ→ℓ+1 (ûℓ ) . All level networks {θℓ }L−1 ℓ=0 are then optimized jointly against the final high-fidelity target with a relative L2 loss, ∥ûL − uL ∥2 Le2e = . ∥uL ∥2 + ϵ This stage does not change the cascade architecture or the one-step inference budget; it adapts the independently trained per-level refiners to the distribution of conditioning states they actually encounter when composed, and turns the one-step rollout into the optimized operating point rather than a coarse approximation to a many-step ODE solve. Section 5.5 confirms empirically that increasing the per-level NFE on the fine-tuned cascade does not improve NRMSE. The calibrated source additionally defines an optional stochastic sampling interface (sample εℓ ∼ qδℓ at each cascade level and integrate the learned velocity field) that we use only as a diagnostic of predictive spread; a stochastic end-to-end variant is reported in Table 7 (Appendix C.4). 7
3.6 Connection to multigrid and inference cost The cascade has a direct analogy to multigrid in classical numerical analysis. In multigrid, a coarse-to-fine V-cycle alternates relaxation steps at multiple grid resolutions to accelerate iterative-solver convergence: the coarsest grid handles low-frequency error modes that fine-grid relaxation cannot reach efficiently, and successive finer grids correct progressively higher-frequency components. MFFM’s cascade is a learned refinement analog: each level corrects residual frequencies that the prolongation from the coarser level fails to capture, and the deterministic e2e operating point replaces the relaxation step at each level with a single network forward pass. MFFM inherits multigrid’s parameter efficiency: per-level networks can be smaller because each handles a narrower band of residual frequencies, and per-level training is the analog of independent local relaxation operators. The total inference cost is L · cnet per query, where cnet is one velocity-network forward pass and L is the cascade depth. Compared with running the HF solver, the speedup is cHF /(L · cnet ). This is favorable when the HF solver is significantly more expensive than the network, which holds for most parametric PDE workflows of practical interest: HF cost grows with grid resolution and physical complexity, often by orders of magnitude, while network cost is fixed by the architecture and grows slowly with input resolution. Architecture (details in Appendix B). The velocity network is a U-Net-style residual convolutional network that takes the channel concatenation (δt , uLF ) as input and outputs a velocity field on the same lattice; the cascade prolongation Iℓ→ℓ+1 is bilinear. Architecture is not a contribution of this work and operator backbones (FNO, DeepONet, neural-operator transformers) can be substituted.
4. Related work Multi-fidelity surrogates. The classical multi-fidelity decomposition that MFFM adopts traces to Kennedy and O’Hagan (2000) and its nonlinear extensions (Perdikaris et al., 2017), surveyed in Peherstorfer et al. (2018). Classical Gaussian-process variants retain a posterior by construction but scale poorly to field-valued solutions; neural variants such as multifidelity DeepONets (Lu et al., 2022; Howard et al., 2023), infinite-fidelity coregionalization (Li et al., 2022b), and multi-resolution operator architectures (Li et al., 2024, 2022a, 2020) model the LF-to-HF correction as a deterministic regressor. Message-passing PDE solvers (Brandstetter et al., 2022) pursue the deterministic-operator route on graph-structured discretizations. Generative refinement for PDE solutions. Generative approaches to PDE problems (Huang et al., 2024; Bastek et al., 2025; Baldan et al., 2026; Hou et al., 2026) typically draw the source from uninformative noise. Closer super-resolution and forecasting precedents include physics-informed score-matching for inverse problems (Holzschuh et al., 2023), dynamics-informed diffusion forecasting (Rühling Cachay et al., 2023), and high-fidelity flow reconstruction (Shu et al., 2023), each of which uses a noise source above a low-resolution conditioning observation. The closest conceptual neighbor PDE-Refiner (Lippe et al., 2023) inserts a diffusion-based refinement on top of an autoregressive prediction with the coarse predictor as a fixed denoising start, but does not calibrate the source to empirical LF-to-HF residual statistics. Cascade structure and source design. The cascade structure draws on cascaded diffusion (Ho et al., 2022; Saharia et al., 2022) and the deterministic-regime analogs 8
of multi-level Monte Carlo (Giles, 2015) and learned multigrid (Greenfeld et al., 2019), but after deterministic end-to-end fine-tuning MFFM operates each level with a single forward pass. Among flow-matching variants that admit non-trivial source-target couplings, including multisample flow matching (Pooladian et al., 2023), bridge methods (De Bortoli et al., 2021; Shi et al., 2023), and stochastic interpolants (Albergo et al., 2025), none treats the source statistics (residual scale and local correlation) as data-estimable design parameters calibrated to a specific application; MFFM identifies multi-fidelity refinement as the regime where this calibration pays off, since the LF-HF data correlation directly determines a residual source design that makes deterministic cascade refinement effective.
5. Experiments 5.1 Setup Benchmarks. We evaluate on eight benchmarks. Two are super-resolution problems: 2D Darcy flow with log-Gaussian permeability, and 1D Burgers’ equation as a 2D space–time field with both axes in the resolution hierarchy. The remaining six are spatiotemporal forecasting tasks: Shallow Water (SW) and Diffusion Reaction (DR) from PDEBench (Takamoto et al., 2022); Shear Flow tracer (Shear-T), Shear Flow pressure (Shear-P), and Active Matter concentration (AM) from The Well (Ohana et al., 2024); and 2D Navier–Stokes (NS) from the FNO dataset (Li et al., 2021). For the forecasting tasks the input is a block of T = 25 past frames at the coarse grid and the target is the corresponding block of T = 25 future frames at the high-fidelity grid, with the temporal axis stacked as channels. Methods compared. We compare nine methods. Bilinear performs no learning and upsamples uLF to the HF grid (Section 5.4 discusses its role as a task-difficulty zero line). FNO-direct (Li et al., 2021) and DeepONet-direct (Lu et al., 2021) take the PDE parameter ξ as input and predict uHF directly. FNO-residual and DeepONet-residual take uLF as input and predict uHF with the same conditioning as MFFM but a deterministic operator backbone. F-FNO (Tran et al., 2023) and CNO (Raonic et al., 2023) are recent neural operators in the same architectural family. PDE-Refiner (Lippe et al., 2023) is the closest probabilistic baseline (iterative diffusion-style refinement; K = 4 steps). FM (single-level) is the single-shot variant of MFFM with L = 1 at the largest parameter count we tested, isolating the cascade contribution. MFFM (cascade) is the method of Section 3. Metrics and protocol. We report global NRMSE; all MFFM main numbers use the deterministic end-to-end cascade point predictor. All methods are trained with AdamW at learning rate 10−5 , cosine decay, and gradient clipping at norm 1.0. MFFM additionally uses EMA decay 0.999 during per-level pretraining and an end-to-end fine-tuning stage at one fifth of the base learning rate, with checkpoint selection by validation NRMSE. 5.2 Main results Table 1 reports NRMSE on all eight benchmarks for the nine methods, split into the two super-resolution problems and the six spatiotemporal forecasting tasks. MFFM-cascade achieves the best NRMSE on seven of the eight benchmarks. The exception is Darcy, where F-FNO obtains the lowest error (0.05830 vs. MFFM-cascade 9
Table 1: NRMSE on eight PDE benchmarks (mean ± std over 3 seeds). For readability the table is split into two halves: top, Darcy, Burgers, SW, DR; bottom, Shear-T, Shear-P, AM, NS. Darcy and Burgers are spatial super-resolution problems; the other six are spatiotemporal forecasting (T =25 past LF frames → future HF frames). Bold marks the best method per benchmark. FNO/DeepONet (direct) take ξ as input; all other learning methods take the LF solution. Method
Darcy
Bilinear 0.08616 FNO (direct) 0.1994 ± 0.0041 FNO (residual) 0.07051 ± 0.024 DeepONet (direct) 0.2467 ± 0.011 DeepONet (residual) 0.06602 ± 0.00031 F-FNO 0.05830 ± 0.00040 CNO 0.8090 ± 0.42 PDE-Refiner 0.09255 ± 0.014 MFFM (cascade)
Method
0.06198 ± 0.0011
Burgers
SW
DR
0.1541 0.1099 ± 0.022 0.06595 ± 0.0057 0.6291 ± 0.0014 0.1269 ± 0.0010 0.08879 ± 0.016 0.1376 ± 0.028 0.1439 ± 0.0040
0.1007 0.05037 ± 0.0027 0.03905 ± 0.0033 0.08304 ± 0.000042 0.1058 ± 0.0000040 0.03575 ± 0.00093 0.1596 ± 0.020 0.1351 ± 0.00029
3.170 0.5557 ± 0.0047 0.7910 ± 0.012 0.9655 ± 0.0034 1.8302 ± 0.0019 1.102 ± 0.18 3.215 ± 0.25 2.378 ± 0.020
0.02488 ± 0.00049 0.004756 ± 0.000076 0.2509 ± 0.0042
Shear-T
Shear-P
AM
NS
Bilinear FNO (direct) FNO (residual) DeepONet (direct) DeepONet (residual) F-FNO CNO PDE-Refiner
0.6260 0.3333 ± 0.0046 0.3278 ± 0.0037 0.7080 ± 0.00021 0.5066 ± 0.0018 0.3161 ± 0.0074 0.3928 ± 0.0023 0.5854 ± 0.0016
0.9011 0.9372 ± 0.0045 3.566 ± 0.18 0.9799 ± 0.011 8.491 ± 0.0035 2.550 ± 0.44 13.89 ± 1.7 59.50 ± 3.2
0.003635 0.9672 ± 0.012 0.02644 ± 0.0051 0.007937 ± 0.00063 0.003441 ± 0.000012 0.01311 ± 0.010 0.09396 ± 0.015 0.09890 ± 0.00030
0.6172 0.1919 ± 0.0040 0.1952 ± 0.0041 0.6104 ± 0.013 0.4185 ± 0.0094 0.1984 ± 0.0067 0.2404 ± 0.010 0.6125 ± 0.000055
MFFM (cascade)
0.2253 ± 0.0012 0.2671 ± 0.0033 0.003162 ± 0.000017 0.04779 ± 0.00031
0.06198); this is also the benchmark where bilinear is already strong (0.08616), leaving less headroom for refinement. Among the deterministic operator baselines, FNO and F-FNO are the strongest competitors but inconsistent across datasets: FNO-direct is competitive on Diffusion Reaction but fails on Active Matter, FNO-residual degrades sharply on Shear-P, and F-FNO wins Darcy but stays well behind MFFM-cascade on the harder tasks. The catastrophic-divergence pattern is most visible on Shear-P (bilinear NRMSE > 0.9): FNO-residual (3.566), CNO (13.89), and PDE-Refiner (59.50) all fail to fit, while MFFM-cascade obtains 0.2671. The largest absolute improvements occur where the LF-to-HF gap is large. On Shallow Water, MFFM-cascade reaches 0.004756 vs. 0.03905 for FNO-residual and 0.04761 for single-level FM (Table 2); on Diffusion Reaction, 0.2509 vs. 0.7910 for FNO-residual; on Navier–Stokes, 0.04779 vs. ∼ 0.20 for the FNO variants and 0.3136 for single-level FM (Figure 2, Appendix C.5, shows a qualitative example). The single-level ablation in Table 2 isolates the cascade contribution: MFFM-cascade improves on single-level FM on every benchmark and by an order of magnitude on Burgers, SW, DR, and NS, supporting the claim that the adjacent-fidelity decomposition is structural rather than a side-effect of the residual flow-matching objective. 10
Table 2: Ablation of MFFM components (global NRMSE). “w/o Gaussian blur”: diagonal source without local correlation. “IID noise”: isotropic noise of matched average variance. “MFFM-single”: one residual flow from coarsest to finest grid (no cascade). “FM-field”: flow matching on the high-fidelity field rather than the residual. Method
Darcy
Burgers
SW
DR
Shear-T Shear-P
MFFM-cascade 0.06198 0.02488 0.004756 0.2509 0.2253 w/o Gaussian blur 0.07562 0.03591 0.004633 0.2701 0.2363 IID noise 0.1365 0.0396 0.02062 0.5081 0.4008 MFFM-single 0.08182 0.1091 0.04761 0.6871 0.5091 FM-field 0.5360 1.026 0.02062 0.5828 0.4848
AM
NS
0.2671 0.003162 0.04779 0.4358 0.003173 0.05224 0.7120 0.003214 0.1151 0.7057 0.003819 0.3136 0.5996 0.9578 0.3051
5.3 Ablation Studies Table 2 ablates the four main design choices of MFFM: the calibrated residual source (vs. IID noise of matched average variance), local Gaussian-blur correlation (vs. a diagonal source), the multi-resolution cascade (vs. a single-shot model), and residual flow matching (vs. flow matching on the high-fidelity field). Each component contributes. Replacing the calibrated source with IID noise causes the largest and most consistent source-side degradation (e.g., 0.04779 → 0.1151 on Navier–Stokes, 0.2671 → 0.7120 on Shear-P), indicating that matching residual statistics matters more than matching average noise scale; Gaussian-blur correlation is a smaller but consistent additional gain on seven of eight benchmarks. Removing the cascade (MFFM-single) produces the largest drop on benchmarks with large LF-to-HF gaps (Burgers, SW, DR, NS). FM-field is substantially worse than residual MFFM on most benchmarks, confirming that flow matching is most effective here when applied to the unresolved bilinear residual rather than to the full high-fidelity field. A control variant that exposes the single-level model to multi-resolution training pairs (Table 6, Appendix C.3) does not close the gap to the full cascade, so the cascade gain is compositional rather than a side-effect of multi-resolution data exposure. Ranking the four design choices by their effect on accuracy across the eight benchmarks: residual decomposition vs. field-space flow matching produces the largest swing, by an order of magnitude on Darcy, Burgers, and Active Matter; the cascade structure is the next largest contributor on the harder benchmarks; calibrated source statistics dominate on the easier benchmarks where the cascade gain saturates; Gaussian-blur correlation is the smallest of the four. 5.4 Bilinear as the task-difficulty zero line The bilinear baseline performs no learning: it upsamples uLF to the HF grid and reports the resulting NRMSE, which equals the relative L2 norm of the residual ∥δ∥ / ∥uHF ∥. Reading the bilinear rows of Table 1 as a task-difficulty axis organizes the benchmarks into a low-residual regime (Active Matter, Darcy, Shallow Water), an intermediate-residual regime (Burgers, Shear-T, Navier–Stokes), and a large-residual regime (Diffusion Reaction, Shear-P) where the LF solution is qualitatively wrong at the HF grid. Across this axis, MFFM remains competitive when bilinear is already strong, improves sharply when the LF-HF gap carries 11
learnable local structure, and recovers useful predictions in the regime where several operator baselines degrade. Per-regime numerical breakdown appears in Appendix C.1. 5.5 Inference cost The deterministic MFFM-cascade used in Table 1 performs one velocity-network evaluation per cascade level. This is the inference mode optimized by deterministic end-to-end finetuning: at each level the residual state is initialized at zero and one midpoint evaluation is applied. Thus, a cascade with L refinement levels requires exactly L velocity-network evaluations for point prediction. Table 5 (Appendix C.2) reports a representative NFE-scaling study in which trained checkpoints are evaluated with NFE/level ∈ {1, 2, 5, 10, 50}. Increasing NFE does not improve NRMSE and only increases cost, because multi-step integration evaluates the velocity networks along off-training deterministic trajectories. We therefore use NFE/level = 1 for all main deterministic results.
6. Conclusion We presented Multi-Fidelity Flow Matching, a refinement framework for high-fidelity PDE solutions that conditions on a low-fidelity solve and uses a residual-calibrated source for flow-matching pretraining. Conditioning makes the LF-to-HF residual substantially easier to model than the full HF field, and the adapted source matches the scale and local structure of empirical PDE residuals. After independent level-wise pretraining, deterministic end-to-end cascade fine-tuning turns the model into a one-step-per-level refiner that reaches the finest grid in L network evaluations per query. Three caveats bound the scope. (i) MFFM requires paired LF and HF training data. (ii) The one-step operating point is optimized for deterministic refinement; stochastic dynamics or weaker LF-HF correlation may require richer source distributions or stochastic end-toend training. (iii) Pretraining and inference use different conditioning chains; end-to-end fine-tuning addresses this empirically, and a formal stability guarantee is future work.
12
References Michael S. Albergo, Nicholas M. Boffi, and Eric Vanden-Eijnden. Stochastic interpolants: A unifying framework for flows and diffusions. Journal of Machine Learning Research (JMLR), 26, 2025. Giacomo Baldan, Qiang Liu, Alberto Guardone, and Nils Thuerey. Physics vs distributions: Pareto optimal flow matching with physics constraints. In International Conference on Learning Representations (ICLR), 2026. Jan-Hendrik Bastek, WaiChing Sun, and Dennis M. Kochmann. Physics-informed diffusion models. In International Conference on Learning Representations (ICLR), 2025. Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural PDE solvers. In International Conference on Learning Representations (ICLR), 2022. Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems (NeurIPS), 2018. Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet. Diffusion Schrödinger bridge with applications to score-based generative modeling. In Advances in Neural Information Processing Systems (NeurIPS), 2021. Alexander I. J. Forrester, András Sóbester, and Andy J. Keane. Multi-fidelity optimization via surrogate modelling. Proceedings of the Royal Society A, 463(2088):3251–3269, 2007. Michael B. Giles. Multilevel Monte Carlo methods. Acta Numerica, 24:259–328, 2015. Daniel Greenfeld, Meirav Galun, Ronen Basri, Irad Yavneh, and Ron Kimmel. Learning to optimize multigrid PDE solvers. In International Conference on Machine Learning (ICML), 2019. Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems (NeurIPS), 2020. Jonathan Ho, Chitwan Saharia, William Chan, David J. Fleet, Mohammad Norouzi, and Tim Salimans. Cascaded diffusion models for high fidelity image generation. Journal of Machine Learning Research (JMLR), 2022. Benjamin Holzschuh, Simona Vegetti, and Nils Thuerey. Solving inverse physics problems with score matching. In Advances in Neural Information Processing Systems (NeurIPS), 2023. Xianglong Hou, Xinquan Huang, and Paris Perdikaris. CFO: Learning continuous-time PDE dynamics via flow-matched neural operators. In International Conference on Learning Representations (ICLR), 2026. Amanda A. Howard, Mauro Perego, George E. Karniadakis, and Panos Stinis. Multifidelity deep operator networks for data-driven and physics-informed problems. Journal of Computational Physics, 2023. 13
Jiahe Huang, Guandao Yang, Zichen Wang, and Jeong Joon Park. Diffusionpde: Generative pde-solving under partial observation. In Advances in Neural Information Processing Systems (NeurIPS), 2024. Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models. In Advances in Neural Information Processing Systems (NeurIPS), 2022. Marc C. Kennedy and Anthony O’Hagan. Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87(1):1–13, 2000. Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research (JMLR), 24(89): 1–97, 2023. Shibo Li, Wei Xing, Robert Kirby, and Shandian Zhe. Multi-fidelity Bayesian optimization via deep neural networks. In Advances in Neural Information Processing Systems (NeurIPS), 2020. Shibo Li, Zheng Wang, Robert Kirby, and Shandian Zhe. Deep multi-fidelity active learning of high-dimensional outputs. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2022a. Shibo Li, Zheng Wang, Robert Kirby, and Shandian Zhe. Infinite-fidelity coregionalization for physical simulation. In Advances in Neural Information Processing Systems (NeurIPS), 2022b. Shibo Li, Xin Yu, Wei Xing, Mike Kirby, Akil Narayan, and Shandian Zhe. Multi-resolution active learning of Fourier neural operators. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2024. Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations (ICLR), 2021. Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. In International Conference on Learning Representations (ICLR), 2023. Phillip Lippe, Bas Veeling, Paris Perdikaris, Richard E. Turner, and Johannes Brandstetter. Pde-refiner: Achieving accurate long rollouts with neural pde solvers. In Advances in Neural Information Processing Systems (NeurIPS), 2023. Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. In International Conference on Learning Representations (ICLR), 2023. 14
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George E. Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Intelligence, 2021. Lu Lu, Raphael Pestourie, Steven G. Johnson, and Giuseppe Romano. Multifidelity deep neural operators for efficient learning of partial differential equations. Physical Review Research, 2022. Ruben Ohana, Michael McCabe, Lucas Meyer, Rudy Morel, Fruzsina J Agocs, Miguel Beneitez, Marsha Berger, Blakesley Burkhart, Stuart B Dalziel, Drummond B Fielding, et al. The well: a large-scale collection of diverse physics simulations for machine learning. Advances in Neural Information Processing Systems, 37:44989–45037, 2024. Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Review, 60(3):550–591, 2018. Paris Perdikaris, Maziar Raissi, Andreas Damianou, Neil D. Lawrence, and George E. Karniadakis. Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling. Proceedings of the Royal Society A, 2017. Aram-Alexandre Pooladian, Heli Ben-Hamu, Carles Domingo-Enrich, Brandon Amos, Yaron Lipman, and Ricky T. Q. Chen. Multisample flow matching: Straightening flows with minibatch couplings. In International Conference on Machine Learning (ICML), 2023. Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci, Rima Alaifari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning of PDEs. In Advances in Neural Information Processing Systems (NeurIPS), 2023. Salva Rühling Cachay, Bo Zhao, Hailey Joren, and Rose Yu. DYffusion: A dynamicsinformed diffusion model for spatiotemporal forecasting. In Advances in Neural Information Processing Systems (NeurIPS), 2023. Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Burcu Karagol Ayan, S. Sara Mahdavi, Rapha Gontijo Lopes, et al. Photorealistic text-to-image diffusion models with deep language understanding. In Advances in Neural Information Processing Systems (NeurIPS), 2022. Yuyang Shi, Valentin De Bortoli, Andrew Campbell, and Arnaud Doucet. Diffusion Schrödinger bridge matching. In Advances in Neural Information Processing Systems (NeurIPS), 2023. Dule Shu, Zijie Li, and Amir Barati Farimani. A physics-informed diffusion model for high-fidelity flow field reconstruction. Journal of Computational Physics, 2023. Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In International Conference on Learning Representations (ICLR), 2021. 15
Yifei Sun, Tao Wang, Junqi Qu, Yushun Dong, Hewei Tang, and Shibo Li. COMPOL: A unified neural operator framework for scalable multi-physics simulations. arXiv preprint, 2025. Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. PDEBench: An extensive benchmark for scientific machine learning. In Advances in Neural Information Processing Systems (NeurIPS) Datasets and Benchmarks Track, 2022. Alexander Tong, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector-Brooks, Kilian Fatras, Guy Wolf, and Yoshua Bengio. Improving and generalizing flow-based generative models with minibatch optimal transport. Transactions on Machine Learning Research (TMLR), 2024. Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized Fourier neural operators. In International Conference on Learning Representations (ICLR), 2023.
16
Appendix A. Conditioning and source-scale analysis This appendix expands the mechanism discussed in Section 3.3. We separate two roles that are coupled in MFFM: conditioning on the low-fidelity solution, which makes the residual refinement problem structurally easier, and residual-calibrated source design, which normalizes the flow-matching regression target. The main experiments use a deterministic endto-end fine-tuned cascade for point prediction; the analysis below explains why the residual flow-matching pretraining provides a favorable initialization for that one-step operating point. A.1 Conditioning, not matched marginals, is the central mechanism The conditional flow-matching loss regresses the velocity network onto the per-sample conditional velocity δ − ε evaluated at interpolated states δt = tδ + (1 − t)ε. The Bayes-optimal velocity represented by the trained network is the conditional regression v ⋆ (z, t, uLF ) = E[δ − ε | δt = z, t, uLF ]. The key structural simplification comes from conditioning on uLF . In the deterministic multi-fidelity setting considered here, the high-fidelity solution is a deterministic function of the underlying PDE instance. If the low-fidelity solution identifies that instance sufficiently well, then the conditional law of the residual collapses around a deterministic refinement target g(uLF ) := uHF − uLF . In the idealized injective case, this conditional law is a Dirac measure at g(uLF ). Given uLF and t < 1, the source can be recovered from a point on the interpolation path by ε=
δt − t g(uLF ) . 1−t
Substituting this identity into the target δ − ε gives the conditional Bayes regressor v ⋆ (z, t, uLF ) =
g(uLF ) − z . 1−t
Along the stochastic interpolation path initialized at δ0 = ε, this velocity reduces to the constant direction g(uLF ) − ε. This idealized calculation explains why a strongly conditioned residual flow is much easier than unconditional field generation. The unconditional marginal velocity does not have this property. For example, if δ and ε are independent zero-mean Gaussians with equal covariance, then ⋆ vmarginal (z, t) = E[δ − ε | δt = z, t] =
2t − 1 z, 2t2 − 2t + 1
(9)
which is explicitly time-dependent. The conditional and marginal calculations are not in conflict: the former exploits the information in uLF , whereas the latter does not. In the main experiments, the final point predictor is not the stochastic path above. After level-wise flow-matching pretraining, we fine-tune the composed cascade using a deterministic 17
2 Table 3: Per-coordinate residual variance ratio σδ2 /σHF from Eq. (11), as a function of empirical correlation ρ between fidelities.
ρ
0.50
0.70
0.90
0.95
0.99
2 σδ2 /σHF
1.00
0.60
0.20
0.10
0.02
rollout with zero residual initialization and one midpoint velocity evaluation per level. Thus, the calculation above should be read as a conditioning argument for why residual flow matching is a good pretraining objective, not as a claim that the final deterministic e2e predictor is a classical many-step ODE solver. A.2 Residual-calibrated sources as normalization The source distribution in MFFM is a design choice. The diagonal residual source 2 d D [δj ], εdiag ∼ N 0, diag(σδ2 ) , σδ,j = Var matches the empirical per-coordinate residual scale. For analytic clarity, consider this diagonal case. If εdiag is independent of δ and has the same per-coordinate variance as the empirical residual, then the regression target satisfies E∥δ − εdiag ∥22 = 2 Tr diag(σδ2 ) . This avoids the Θ(n) target scale induced by an isotropic unit-variance source on unitnormalized fields and keeps the flow-matching regression problem numerically well scaled. The main implementation further applies local Gaussian-blur correlation. Concretely, we sample iid noise, apply a depthwise Gaussian blur, normalize each sample to unit empirical standard deviation, and then scale by σδ . This preserves the residual scale used in the diagonal analysis while adding local spatial correlation, which better matches the structure of PDE residual fields. The ablations in Section 5.3 show that removing Gaussian blur or replacing the residual-calibrated source with IID noise degrades accuracy. A.3 Per-coordinate variance reduction Let ρj = Corr(uLF [j], uHF [j]) denote the empirical correlation between the two fidelities at coordinate xj . The residual variance is 2 2 2 σδ,j = Var[uHF [j] − uLF [j]] = σHF,j + σLF,j − 2ρj σHF,j σLF,j .
(10)
When the two fidelities have similar per-coordinate variance (σLF,j ≈ σHF,j ), this simplifies to 2 2 σδ,j ≈ 2(1 − ρj )σHF,j . (11) Table 3 reports the corresponding variance ratio for representative correlations. This calculation explains why adjacent-fidelity residual refinement is easier than full-field generation when fidelities are strongly correlated. It also motivates the cascade: adjacent levels typically have higher correlation than the coarsest-to-finest pair, reducing the scale of the residual that each level must learn. 18
Table 4: Per-benchmark data and training configuration. “Split” reports train/validation/test sample counts. “Pretrain epochs” denotes independent level-wise flow-matching training. “E2E FT epochs” denotes deterministic end-to-end cascade fine-tuning after level-wise pretraining. Cascade resolutions list the grids visited by the level-wise refinement. Benchmark Darcy Burgers Shallow Water Diffusion Reaction Shear Flow (tracer) Shear Flow (pressure) Active Matter Navier–Stokes
Task type
Cascade resolutions
Split
spatial SR space–time SR block forecasting block forecasting block forecasting block forecasting block forecasting block forecasting
{16, 32, 64, 128}2
512/128/128 512/128/128 800/100/100 800/100/100 896/112/112 896/112/112 175/24/26 2048/256/256
{32, 64, 128, 256}2 {16, 32, 64, 128}2 {16, 32, 64, 128}2 {16, 32, 64, 128}2 {16, 32, 64, 128}2 {16, 32, 64, 128}2 {8, 16, 32, 64}2
Pretrain epochs E2E FT epochs 500 100 100 300 300 300 100 300
100 100 100 300 300 300 100 100
A.4 Relation to deterministic one-step inference The analysis above does not imply that the final e2e fine-tuned model should be evaluated with many ODE steps. In fact, the main deterministic predictor is explicitly trained as a one-step-per-level cascade. At level ℓ, the rollout used during e2e fine-tuning is 1 ûℓ+1 = Iℓ→ℓ+1 (ûℓ ) + vθℓ 0, , Iℓ→ℓ+1 (ûℓ ) . 2 Therefore, increasing the number of Euler evaluations at test time changes the operating distribution of the velocity networks. The deterministic NFE-scaling study in Appendix C.2 confirms that larger NFE does not improve the fine-tuned predictor. This is why MFFMcascade is used as a one-step-per-level deterministic refiner in all main results.
Appendix B. Implementation details This appendix collects the architectural, data-construction, and training-protocol details elided from Section 3.6. B.1 Per-benchmark configuration Table 4 lists the cascade resolutions, train/validation/test split sizes, and training-stage epoch budgets used for each benchmark. The same shared training protocol described in Section 5.1 (AdamW, cosine schedule, gradient clipping, EMA) is used throughout; the only per-benchmark variation is the resolution hierarchy, the split sizes, and the epoch budgets reported here. B.2 Data construction We use two data-construction protocols: super-resolution of paired solution fields and spatiotemporal block forecasting. In all cases, lower-resolution fields are bilinearly prolonged to the target resolution before forming residuals. Darcy. The Darcy benchmark is generated in-house as a multi-resolution elliptic PDE dataset. We sample 768 coefficient fields from a Gaussian random field prior and construct 19
lognormal permeability fields. For each sample, the same underlying high-resolution coefficient field is downsampled to the cascade resolutions {16, 32, 64, 128}2 , and the Darcy equation is solved separately on each grid. Thus, u16 , u32 , u64 , and u128 are paired multi-resolution solver outputs for the same underlying PDE instance, rather than resized copies of the highresolution solution. The resulting paired fields are split into 512/128/128 train/validation/test samples. The cascade learns spatial residual refinements over the hierarchy {16, 32, 64, 128}2 . Burgers. The Burgers benchmark is generated in-house from the one-dimensional viscous Burgers equation. Although the underlying PDE is one-dimensional in space, we represent each solution trajectory as a two-dimensional space–time field u(t, x). For each of 768 samples, we draw one continuous Gaussian-random-field initial condition and solve the same underlying sample at resolutions {32, 64, 128, 256}. At resolution s, we use s spatial points and save s future time snapshots, giving an output tensor of shape s × s. Thus, both the spatial and temporal axes are included in the resolution hierarchy, and the task is space–time super-resolution rather than one-step temporal forecasting. The data are split into 512/128/128 train/validation/test samples. Spatiotemporal forecasting benchmarks. For Shallow Water, Diffusion Reaction, Shear Flow, Active Matter, and Navier–Stokes, we use block-to-block forecasting. Each sample contains 50 frames. The first 25 frames are used as the coarse-grid input block, and the last 25 frames are used as the high-fidelity target block. The temporal axis is stacked into the channel dimension, so a block with 25 frames and Cphys physical variables is represented as a tensor with 25Cphys channels. Shallow Water and Diffusion Reaction are taken from PDEBench; Shear Flow and Active Matter are taken from The Well; Navier–Stokes is taken from the FNO dataset PDENavierStokes_V1e-3_N5000_T50. When the native spatial grid is not square, we resize the fields to square grids before constructing the cascade resolutions. The level-wise residual at cascade level ℓ is written in a form that covers both super-resolution and past-to-future block forecasting: δℓ = utar ℓ+1 − Iℓ→ℓ+1 (cℓ ), where utar ℓ+1 denotes the target field or target frame block at resolution level ℓ + 1, and cℓ denotes the conditioning field at level ℓ. For ℓ = 0, c0 is the observed coarse input field or coarse past-frame block. For later cascade levels, cℓ is the lower-resolution target field during level-wise pretraining and the prediction propagated from the previous cascade level during inference. For Darcy and Burgers, the conditioning field is simply the lower-resolution solution field, cℓ = uℓ , so the residual reduces to δℓ = uℓ+1 − Iℓ→ℓ+1 (uℓ ). For the spatiotemporal forecasting benchmarks, utar ℓ+1 is the stacked future-frame target block, while cℓ is the corresponding conditioning block used at that cascade level. B.3 Architecture rationale The velocity network is a U-Net-style residual convolutional network mapping the channel concatenation (δt , uLF ) ∈ RB×2C×H×W to a velocity field in RB×C×H×W on the same lattice. 20
Here C denotes the number of solution channels, equal to T · Cphys in the spatiotemporal forecasting setting where time frames are stacked as channels. The network uses local convolutions, residual blocks, skip additions, and time-dependent FiLM conditioning. We use this lightweight image-to-image backbone because the architectural choice is not the main contribution; the method only requires a velocity network that takes (δt , uLF , t) and outputs a residual velocity on the target grid. The architectural extension specific to MFFM is the channel concatenation of uLF with δt at the network input, which exposes the cheap solution to all convolutional layers. We use channel concatenation rather than a global conditioning bottleneck because uLF is a structured spatial field rather than a low-dimensional descriptor. The cascade prolongation Iℓ→ℓ+1 is bilinear interpolation; learned prolongations would couple cascade levels and break the simple level-wise construction. B.4 Velocity network Each per-level velocity network vθ has the following components. Input and output. The input is the channel concatenation (δt , uLF ) ∈ RB×2C×H×W , where δt is the interpolated residual state and uLF denotes the prolongated conditioning field at the current level. The output is a velocity field vθ (δt , t, uLF ) ∈ RB×C×H×W . Input lift. A 3×3 convolution maps the 2C input channels to a hidden width h. Time conditioning. The scalar time t ∈ [0, 1] is embedded with a sinusoidal embedding and passed through a two-layer MLP to produce a hidden time-conditioning vector. Each residual block receives this conditioning vector through a FiLM module. Residual blocks and skip additions. The network applies a stack of residual convolutional blocks followed by a second stack with skip additions from the first stack. Each block uses group normalization, SiLU activations, two 3×3 convolutions, and FiLM modulation from the time embedding. All convolutions operate at the current target-grid resolution; the velocity network itself does not perform spatial downsampling or upsampling. Output head. A final group normalization, SiLU activation, and zero-initialized 1×1 convolution map the hidden features back to C channels. Zero initialization makes the initial velocity approximately zero, giving a stable identity-like initialization for the ODE update. Group normalization. We use group normalization with up to eight groups, choosing a divisor of the hidden width when necessary. Exponential moving average. During independent level-wise flow-matching pretraining, model weights are tracked by an EMA copy with decay 0.999, and the EMA weights are applied before the cascade fine-tuning stage. The deterministic end-to-end fine-tuning stage does not use EMA; the final checkpoint is selected by validation NRMSE. B.5 Training protocol Optimizer. Independent level-wise flow-matching pretraining uses AdamW with learning rate 10−5 , β1 = 0.9, β2 = 0.999, weight decay 10−6 , and gradient clipping at L2 norm 21
1.0. EMA weights with decay 0.999 are used during this pretraining stage. After level-wise pretraining, MFFM-cascade is fine-tuned end-to-end using the deterministic cascade rollout with one fifth of the base learning rate, the same weight decay, and the same gradient clipping. The end-to-end fine-tuned checkpoint is selected by validation NRMSE. Schedule. Cosine learning-rate decay over the full training horizon, no warmup at small scale. Batch and epoch counts. Batch size is 32. Per-benchmark settings appear in Table 4. Source-statistics precomputation. Before training each level, the empirical per-coordinate d D [δj ] is computed on the corresponding training residuals and residual variance σj2 = Var stored. For cascade level ℓ, the residual is δℓ = utar ℓ+1 − Iℓ→ℓ+1 (cℓ ). For pure super-resolution tasks, this reduces to δℓ = uℓ+1 − Iℓ→ℓ+1 (uℓ ). The source noise is generated by sampling iid Gaussian noise, applying a depthwise Gaussian blur kernel, normalizing each sample to unit empirical standard deviation, and scaling by the stored residual standard deviation σδ . The variance is floored at 10−8 to avoid numerical issues at coordinates with near-zero residual variance. Per-level width and depth. The hidden width and number of residual blocks are selected per cascade level. Earlier levels usually handle larger LF-to-HF corrections, while later levels mainly perform smaller high-resolution refinements, so later levels are assigned no larger capacity than earlier ones. In the implementation used for the main experiments, the hidden width is lower-bounded by 32 and rounded to a multiple of eight, and the number of residual blocks is lower-bounded by two. Exact per-level values are recorded in the released training scripts. B.6 Algorithm pseudocode Algorithm 2 gives the deterministic rollout used both for point prediction and for end-to-end cascade fine-tuning. During fine-tuning, the returned ûL is compared with the finest-resolution target, and all level networks are optimized jointly.
Appendix C. Additional empirical analysis C.1 Bilinear-axis regime breakdown This appendix expands the per-regime discussion summarized in Section 5.4. The bilinear baseline measures the relative size of the LF–HF residual before learning, NRMSEbilinear =
∥uHF − I(uLF )∥2 . ∥uHF ∥2
It therefore provides a simple axis for interpreting where refinement methods have room to improve. 22
Algorithm 1 MFFM training (single cascade level) (i)
tar,(i)
Require: Level-wise paired data Dℓ = {(ũℓ , uℓ+1 )}N i=1 , network vθ , total steps T tar,(i)
1: δ (i) ← uℓ+1
(i)
− ũℓ for all training pairs d D [δ] ▷ per-coordinate residual variance 2: σδ2 ← Var ℓ 3: for step = 1 to T do 4: Sample minibatch (ũℓ , utar ℓ+1 ) ∼ Dℓ tar 5: δ ← uℓ+1 − ũℓ 6: Sample z ∼ N (0, I) 7: z̄ ← Kτ ∗ z ▷ depthwise Gaussian blur d 8: z̃ ← z̄/(Std(z̄) + ϵnum ) 9: ε ← σδ ⊙ z̃ 10: Sample t ∼ U[0, 1] 11: δt ← t δ + (1 − t) ε 12: J ← ∥vθ (δt , t, ũℓ ) − (δ − ε)∥2 13: Update θ by AdamW on J ; update EMA weights 14: end for 15: Apply EMA weights to vθ 16: return vθ Algorithm 2 Deterministic cascade point prediction: produces ûL at the finest grid using one velocity-network evaluation per cascade level. Require: Coarse-grid input u0 , trained per-level networks {vθℓ }L−1 ℓ=0 1: û0 ← u0 2: for ℓ = 0 to L − 1 do 3: ũℓ ← Iℓ→ℓ+1 (ûℓ ) 4: δℓ ← 0 5: vℓ ← vθℓ (δℓ , 1/2, ũℓ ) 6: ûℓ+1 ← ũℓ + vℓ 7: end for 8: return ûL
Low-residual regime. Active Matter (0.0036), Darcy (0.0862), and Shallow Water (0.1007) have relatively small bilinear residuals. In this regime, the LF input already captures much of the HF structure, although Shallow Water still contains learnable high-frequency error: MFFM reduces its NRMSE from 0.1007 to 0.004756. On Active Matter and Darcy, the absolute headroom is smaller, and the difference between strong refinement methods is correspondingly modest. Intermediate-residual regime. Burgers (0.1541), Shear Flow tracer (0.6260), and Navier– Stokes (0.6172) require nontrivial local correction beyond interpolation. MFFM improves these to 0.02488, 0.2253, and 0.04779, respectively, indicating that the cascade can recover unresolved local structure from the LF input. 23
Large-residual regime. Diffusion Reaction (3.170) and Shear Flow pressure (0.9011) have the largest bilinear residuals. These are the settings where direct residual learning and operator baselines are most unstable. MFFM reduces the errors to 0.2509 and 0.2671, respectively, supporting the value of splitting the LF–HF correction into adjacent-fidelity residual refinements. C.2 Deterministic NFE scaling The main MFFM-cascade results use the deterministic one-step rollout optimized during endto-end fine-tuning. This appendix tests whether increasing the number of Euler evaluations per cascade level improves the fine-tuned deterministic predictor. For each benchmark, we train a representative single-seed checkpoint with NFE/level = 1 during deterministic end-to-end fine-tuning, then keep all network weights fixed and vary only NFE/level ∈ {1, 2, 5, 10, 50} at evaluation time. The residual state is initialized at zero at each level, matching the deterministic point-prediction mode used in the main experiments. A larger NFE therefore does not correspond to additional training or a different model; it only changes the numerical rollout used at inference.
Table 5: Deterministic NFE-scaling study for MFFM-cascade. Each row uses one representative single-seed checkpoint trained with NFE/level = 1 during deterministic end-to-end fine-tuning; only the number of Euler evaluations per cascade level is changed at evaluation time. The one-step setting is therefore the optimized operating point. Total NFE is NFE/level × L, with L = 3 refinement levels for the listed benchmarks. Table 1 reports multi-seed averages for the main comparison. Benchmark / NFE per level Total NFE Burgers Shallow Water Active Matter Navier–Stokes
1
2
5
10
50
3 0.02556 0.003378 0.003082 0.05996
6 0.03794 0.008787 0.003082 0.1054
15 0.04660 0.01195 0.003084 0.1363
30 0.04900 0.01297 0.003085 0.1458
150 0.05081 0.01376 0.003085 0.1531
Across the tested benchmarks, increasing NFE does not improve the fine-tuned deterministic predictor. The best or tied-best NRMSE is obtained at NFE/level = 1, and larger values either leave the error nearly unchanged or degrade it. This behavior is consistent with the training procedure: end-to-end fine-tuning optimizes the composed one-step rollout directly, while larger NFE evaluates the velocity networks along deterministic trajectories that were not optimized by the fine-tuning objective. The result should therefore be interpreted as an operating-point study rather than as a classical ODE solver convergence test. MFFM-cascade is designed to be used as a onestep-per-level deterministic refiner after end-to-end fine-tuning. Additional Euler evaluations increase the number of velocity-network calls linearly without improving the point predictor. 24
C.3 Multi-resolution training pairs are not the source of the cascade gain Section 5.3 compares the full cascade against a single-level model trained from coarsest to finest. A natural concern is that the cascade improvement reflects exposure to multi-resolution training pairs rather than a compositional benefit of adjacent-fidelity refinement. We control for this with MFFM-single (multires train): a single residual flow trained on the same multi-resolution pairs available to the cascade, but evaluated as a direct coarsest-to-finest refiner. Table 6: Single-level model with multi-resolution training data does not match the full cascade. Numbers are global NRMSE; MFFM-cascade and MFFM-single repeat values from Table 2 for reference. Method
Darcy
Burgers
SW
DR
Shear-T
Shear-P
AM
NS
MFFM-cascade MFFM-single MFFM-single (multires train)
0.06198 0.08182 0.08136
0.02488 0.1091 0.1227
0.004756 0.04761 0.01702
0.2509 0.6871 0.5591
0.2253 0.5091 0.5892
0.2671 0.7057 0.5023
0.003162 0.003819 0.04354
0.04779 0.3136 0.2624
The multires-trained single-level model improves over the ordinary single-level model on Shallow Water and Diffusion Reaction but remains far worse than the cascade; on Shear-T and Active Matter it is even worse than the ordinary single-level variant. The cascade’s compositional adjacent-resolution refinement is therefore the load-bearing structural choice, not multi-resolution data exposure alone. C.4 Stochastic end-to-end cascade variant The main MFFM-cascade results in Table 1 use deterministic end-to-end fine-tuning and deterministic point prediction. We additionally test a stochastic end-to-end variant to assess whether the source-sampling interface can be optimized directly. This experiment is not used for the main accuracy comparison; it is included as an appendix diagnostic. Starting from the same level-wise residual flow-matching pretraining, we fine-tune the composed cascade with stochastic rollouts. At each cascade level, the residual state is initialized by a source sample δℓ,0 = εℓ , εℓ ∼ qδℓ , and the learned velocity field is integrated with one Euler evaluation per level. For a minibatch, we draw Ktrain independent stochastic rollouts and optimize the Monte Carlo estimate of the expected relative error: Lstoch-e2e =
KX train
1 Ktrain
k=1
(k)
∥ûL − uL ∥2 . ∥uL ∥2 + ϵ
In the implementation used here, Ktrain = 4 and Kval = 8. At test time, we draw Ktest = 10 samples and report the NRMSE of the ensemble mean together with the average predictive spread: Ktest 1 X (k) (k) ūL = ûL , s = meanx Stdk [ûL (x)]. Ktest k=1
25
Table 7: Stochastic end-to-end cascade variant. The model is fine-tuned using stochastic cascade rollouts and evaluated by the ensemble mean over Ktest = 10 samples. “Spread” is the mean pointwise standard deviation across stochastic samples. This experiment is a diagnostic of the stochastic sampling interface and is not used for the main deterministic comparison. Benchmark
NRMSE (NFE/level=1)
Darcy Burgers Shallow Water Diffusion Reaction Shear-P Shear-T Active Matter Navier–Stokes
0.07701 0.03274 0.01207 0.3325 0.3052 0.2524 0.003226 0.05711
Spread (NFE/level=1) −4
2.232×10 2.209×10−3 1.012×10−2 1.269×10−2 4.888×10−3 2.897×10−2 1.768×10−3 4.057×10−2
NRMSE (NFE/level=2)
Spread (NFE/level=2)
0.07693 0.04089 0.02815 0.3776 0.5339 0.4175 0.003228 0.1227
2.233×10−4 2.955×10−3 9.940×10−3 1.182×10−2 4.271×10−3 2.300×10−2 1.768×10−3 3.303×10−2
The stochastic variant confirms that the learned residual flow can support a sampling-based predictor, but it is not the operating point used for the main results. In most benchmarks, the stochastic ensemble mean is less accurate than the deterministic e2e predictor in Table 1, especially when moving from one to two Euler evaluations per level. This mirrors the deterministic NFE-scaling behavior in Appendix C.2: the networks are most effective at the one-step operating point for which the cascade is fine-tuned. We therefore use deterministic e2e prediction for the headline NRMSE comparison and treat stochastic rollouts as an optional diagnostic of predictive spread. C.5 Visualization Figure 2 provides a qualitative example on the Navier–Stokes benchmark. We visualize five predicted future frames from one held-out test trajectory. The rows show the MFFM cascade prediction, the high-fidelity ground truth, and the bilinear interpolation baseline, respectively. The columns correspond to representative future-frame indices. This figure is intended as a qualitative diagnostic rather than a separate quantitative claim; the full numerical comparison is reported in Table 1. The example illustrates the role of learned residual refinement. Bilinear interpolation captures the large-scale flow pattern but produces overly smooth structures and misses sharper local variations. The MFFM cascade prediction more closely matches the highfidelity target across the selected frames, consistent with the lower Navier–Stokes NRMSE reported in Table 1.
26
Figure 2: Qualitative Navier–Stokes visualization on a held-out test trajectory. Columns show five selected future frames. Rows show MFFM prediction, high-fidelity ground truth, and bilinear interpolation, respectively. MFFM preserves the large-scale structure while correcting local errors left by pure interpolation.
27