Physics-Constrained Flow Matching: Sampling Generative Models with Hard Constraints
Authors: Utkarsh Utkarsh
Pengfei Cai
Alan Edelman
Rafael Gomez-Bombarelli
Christopher Vincent Rackauckas
Venue: NeurIPS 2025 poster ()
Keywords: generative models
flow matching
scientific machine learning
hard constraints
conservation laws
Primary area: machine_learning_for_sciences

Abstract: Deep generative models have recently been applied to physical systems governed by partial differential equations (PDEs), offering scalable simulation and uncertainty-aware inference. However, enforcing physical constraints, such as conservation laws (linear and nonlinear) and physical consistencies, remains challenging. Existing methods often rely on soft penalties or architectural biases that fail to guarantee hard constraints. In this work, we propose Physics-Constrained Flow Matching (PCFM), a zero-shot inference framework that enforces arbitrary nonlinear constraints in pretrained flow-based generative models. PCFM continuously guides the sampling process through physics-based corrections applied to intermediate solution states, while remaining aligned with the learned flow and satisfying physical constraints. Empirically, PCFM outperforms both unconstrained and constrained baselines on a range of PDEs, including those with shocks, discontinuities, and sharp features, while ensuring exact constraint satisfaction at the final solution. Our method provides a flexible framework for enforcing hard constraints in both scientific and general-purpose generative models, especially in applications where constraint satisfaction is essential.

Physics-Constrained Flow Matching: Sampling Generative Models with Hard Constraints
Utkarsh Massachusetts Institute of Technology
Pengfei Cai* Massachusetts Institute of Technology
Alan Edelman Massachusetts Institute of Technology
Rafael ComeezBommbarele
Massachusetts Institute of Technology
Christopher Vincent Rackauckas Massachusetts Institute of Technology
Abstract
Deep generative models have recently been applied to physical systems governed by partial differential equations (PDEs), offering scalable simulation and uncertainty- aware inference. However, enforcing physical constraints, such conservation laws (linear and nonlinear) and physical consistencies, remains challenging. Exist- ing methods often rely on soft penalties or architectural biases that fail guarantec hard constraints. In this work, we propose (PCFM), zero-shot inference framework that enforces arbitrary nonlinear con- straints in pretrained flow-based generative models. PCFM continuously guides the sampling process through physics- -asee corrections applied to intermediate solu- tion states, while remaining aligned with the learned flow and satisfying physical constraints. Empirically, PCFM outperforms both unconstrained and constrained baselines range of PDEs, including those with shocks, discontinuities, and sharp features, while ensuring exact constraint satisfaction at the final solution. Our method provides flexible framework for enforcing hard constraints in both scientific and general-purpor generative models, especially in applications where constraint satisfaction is essential.
Flow Matching
Introduction
Deep generative modeling provides powerful and data-efficient framework for learning complex distributions from finite samples. By estimating an unknown data distribution and enabling sample generation via latent-variable models, deep generative methods have achieved state-of-the-art per- formance in wide range of domains, including image synthesis (1/4), natural language generation 5 6, and applications in molecular modeling and materials simulation 7-9..
Inspired by these successes, researchers have begun applying generative modeling to physical systems governed by Partial Differential Equations (PDEs) (10H3). In these settings, generative models offer unique advantages, including efficient sampling, uncertainty quantification, and the capacity model multimodal solution distributions [10 12]. However, fundamental challenge this context is ensuring that generated samples respect the governing physical constraints of the system 12 14. In traditional domains like vision or text, domain structure is often incorporated
'Equal contribution. Order decided coin toss. Corresponding authors: raf rafagttmtt edu, crackauc@mit. edu The code implementations are available in Python and Julia
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

through soft constraints clastrinnts-classiffee guidance [3 score conditioning [2], or architectural priors such as equivariance [15]. Manifold-based approaches further constrain generations lie on known geometric spaces [16 While such methods can align the model with geometric priors, they cannot be easily adapted for enforcing physical laws in dynamical systems.
Crucially, constraint enforcement in generative modeling for PDEs follows a different paradigm. Physical invariants such as mass, momentum, and energy [19 20) often arise from underlying symmetries [21]. Prior work to incorporate physics into neural networks has largely relied on
inductive biases in training and regression- based tasks: encoding conservation laws as soft penalties (e.g.,
Neural Networks (PINNs)) (22), imposing architectural priors [23] (2).. However, soft constraints can lead critical failure modes, particularly when exact constraint satisfaction is essential for stability or physical plausibility [25 To address this, recent efforts
have explored hard constraint enforcement, through learning conservation laws 28
constraint
satisfaction at inference [14 30 and differentiable physics (1333)
Despite this progress, hard constraint enforcement in generative models, particularly for PDEs, remains a nascent area [12 34 Enforcing hard constraints in generative models is particularly challenging due the inherent stochasticity of the sampling process, and the constraints must be satisfied exactly in the final denoised solution but need not be preserved throughout the sampling process. DiffusionPDE [0] and D-Flow 351 propose mundient-based constraint enforcement during sampling. but these methods often require backpropagation through expensive PDE operators and may fail exactly satisfy the target constraints. The ECI framework [12] introduces a novel mixing- based correction process for zero-shot constraint satisfaction, but only empirically evaluates on simple linear constraints and shows limited robustness on sharp, high-gradient PDEs with shocks or discontinuities.
In this work, we introduce Physics- Constrained Flow Matching (PCFM), framework that bridges modern generative modeling with classical ideas from numerical PDE solvers. PCFM enables zero-shot constraint enforcement up to machine precision for pretrained flow matching models, projecting intermediate flow states onto constraint man- ifolds at inference time, without requiring gradient information during training. Unlike prior methods, PCFM enforce arbitrary nonlinear equality constraints, including global conservation laws, nonlinear residuals, and sharp boundary conditions. It requires no retraining and no architectural modification, operating entirely post-hoc. While we focus on PDE-constrained generation in this work, PCFM provides a flexible framework for enforcing hard constraints for flow-based generative models and may be extended to work for other domains such molecular design and scientific simulations beyond PDEs. To summarize, we list our contributions follows:
Flowmattchinggsee
Physica- Fiow Matching
Vanilla Flow Matching
Figure 1: Evolution of generated solutions for the Burgers equation using vanilla Flow Match- ing (bottom) and our Physics-Constraimg Matching (top). Burgers' equation exhibits sharp shock fronts (top left in the figure). which standard FFM fails to capture accurately, resulting in overly smoothed smeared solutions. In contrast, PCFM efficiently incorporates physical constraints during sampling, enabling accurate shock resolution and physically consistent final outputs.
Flow
1. We introduce general framework PCFM for enforcing arbitrary and mulriple physical constraints in Flow Matching-based generative models. These constraints include satisfying conservation laws, boundary conditions, even arbitrary non-linear constraints. Our method enforces these constraints as hard requirements at inference time, without modifying the underlying training objective.
2. Our approach is zero-shot: it operates directly on any pre-trained flow matching model without requiring gradient information for the constraints during training. This makes the method broadly applicable and computationally efficient, especially in scenarios where constraint gradients are expensive unavailable.

---

have emerged Table 1: Comparison of generation methods motivated by constraint guidance or enforcement.
Zero-shot
Continuous Guidance
Constraint Satisfaction
Complex Constraints
Conditional FFM 11
X
Dif [10
X
D-Flow 35
X
ECT 12
PCFM (Ours)
3. We demonstrate significant improvements in generating solutions to partial differential equations (PDEs), outperforming state-of-the-art baselines by up to 99.5% in standard metrics, such as mean squared error, while ensuring zero constraint residual.
4. We evaluate our method on challenging PDEs exhibiting shocks, discontinuities, and sharp spikes-setiins in which standard flow matching models typically fail. Our approach improves the accuracy of such models at inference time only, without the need for retraining or fine-tuning, by retrofitting physical consistency into generated samples.
5. To enable practical deployment, we develop custom, batched and differentiable solver that projects intermediate flow states onto the constraint manifold. This solver integrates seamlessly with modern deep learning pipelines and enables end-to-end differentiability through the constraint enforcement mechanism.
2 Related Work
Flow-based Generative Models. Flow-based generative models 44 36 37
scalable alternative to diffusion models by defining continuous normalizing flows (CNF) using ordinary differential equations (ODEs) [38 39 36 parameterized by time-dependent vector field. In flow matching models, samples from a tractable prior distribution are transported to a target distribution via learned vector field through a simulation-free training of CNFs. Stochastic interpolants provide a unifying framework to bridge deterministic flows and stochastic diffusions, enabling flexible interpolations between distributions (37). Furthermore, methods like rectified flows proposed efficient sampling with fewer integration steps by straightening the flow trajectories [36]. Functional flow matching (FFM) [11 and Denoising Diffusion Operator (DDO) 440] extend this paradigm to spatiotemporal data, learning flow models directly over function spaces such as partial differential equations (PDEs) solutions.
Constraint Guided Generation. Guiding generative models with structural constraints is an upcoming direction used to improve fidelity in physics-based based domains [2]
Constraint infor- 41 35 and has been successively applied to domains such as inverse problems 45 Gradient backpropagation through an ODE solver be prohibitively expensive for functional data such PDEs 48) Manifold-based flows and diffusion models [6 [8] capture known geometric priors. However, they are not suitable for PDEs having with data-driven implicit constraints. For PDB-constrained generation, DIF- FUSIONPDE applies PINN-like soft penalties during sampling, while D-FLOw [35] optimizes coide-oonditioned objective. Both approaches incur high computational cost and offer only approximate constraint satisfaction. In the context of scored-based diffusion models, PDM 49
mation is typically exploited via gradient backpropagation [42
10
projects every partially denoised intermediate sample during the sampling phase. However, enforeing constraints at every timestep can over-constrain the trajectory, especially under nonlinear or global introduces
constraints like ICs or mass conservation that only apply the final state. ECI [22 novel gradient-free zero-shot, and hard constraint method on PDE solutions. However, its empirical evaluation is limited to simple linear and non-overlapping constraints e.g., pointwise or regional equalities- wwith known closed-form projections. It lacks a general roadmap for nonlinear or coupled constraints and has limited evaluation on harder-to-solve PDEs with shocks. Furthermore, while labeled "gradient-free", its reliance on analytical projections restricts extensibility to nonlinear or coupled constraints, and practical enforcement still implicitly relies on nonlinear optimization, which often requires gradient information [50]. We summarize our differences compared to generative methods motivated by constraint satisfaction in Table DI
Physics-Informed Learning and Constraint-Aware Numerics. Physics-informed learning in- corporates physical laws as inductive biases in machine learning models, typically through soft penalties as in PINNs 22 5 or Neural Operators [23 52 (33.. While effective for regression, these
methods lack guarantees of hard constraint satisfaction. Hard constraints have been addressed via
3

---

differentiable layers [54] (3] 56] inference- time projections [1, [0], and structured architectures
B2
B3
28
Constraint-aware integration methods offer complementary insights. constrained neural differential algebraic equations [59 60 and geometric integrators [6], 61
ODEs 57
58 29
enforce feasibility through projection- based updates. Though underexplored in generative modeling, these methods motivate principled approaches to constrained sampling. Our method combines this numerical perspective with flow-base bow-based generation, enabling exact constraint enforcement without retraining.
Methodology and Setup
3.1 Problem Setup
We consider physical systems governed by parameterized conservation laws of the form
0,u(x,t)+ V Folu(a, t)) 0, u(x, 0) = ao(2x), Bu(x, t) O,
Vz e Ω t e [0,T], Vx E EΩ, Vx 0n, t E [0,7],
(1)
(2)
(3)
where 52 C Rd is a bounded domain, u x (0, X - R is the solution field, (u)
is flux function parameterized by Φ E Φ, 0g specifies the initial condition, and B denotes the boundary operator. For a fixed PDE family and parameter set 4, we define the associated solution 3 representing all physically admissible solutions generated by varying é We assume access to a pretrained generative model like FFM M.. that approximates this solution set, e.g., via flow-bet www-based model trained on simulated PDE solutions.
set UF == {u Eu: 3Ø E & such that u satisfies [
In addition the governing equations, we consider physical constraint which we wish enforce on the solution u(I, t) through a constraint operator Hu(æ, t) = 0 defined on subdomain N'x S Q x [0, T). and let Uz == uu E U Hu(a, t) for all (xr t) E Xx) denote the constraint- satisfying solution set. Our objective is to generate samples from the intersection UFTH == Uy n nlH, i.e. functions that satisfy both the PDE and the constraint exactly. Importantly, we seek to impose H at inference time, without retraining the pre-trained model, thus narrowing the generative support from Uy to Uyyn in a zero-shot manner.
3.2 Generative Models via Flow-Based Dynamics
Let U denote the space of candidate functions : X - R", where X := 52 [0, is the spatiotemporal domain of the PDE. We define a family of time-dependent diffcomorphic map Φ : u x [o. 1] - 24 called flow, indexed by the flow time T [0. A vector field v/ which defines the evolution of the flow @r according to the ordinary differential equation: o 2, (u))., Φ0(4))= = ug. This yields a continuous path of measures through push-forward map TT (0-) To. connecting the prior noise measure To to the predicted target measure πi. Chen et al. 38 proposed to learn the vector field i2 with a deep neural network parameterized by O. Later in this text, we will denote this parameterized vector field as Uy == vr(x;; f). This naturally induces parameterization for Φ called Continuous Normalizing Flows (CNFs). Lipman et al. [4] introduced a simulation and likelihood free training method for CNFs called Flow Matching along with with other similar works [36] 39
(see Appendix G).
In the context of functional data such as PDE solutions, the Functional Flow Matching (FFM) framework [ extends this idea to functional spaces. Instead of modeling samples in finite- dimensional Euclidean space, FFM learns flows between functions uo, d1 E U via interpolants 4,37 36), trained on trajectories with un To, U1 v, where v is the target probability measure from where wish to sample the PDE solutions. The vector field vo is parameterized by a Neural Operator 23 24 which is forward compatible to operate on functional spaces. The resulting generative model approximates the target measure Ti R V over the support Us through sufficient training (see Appendix G.2 for further details). In this work, we focus on augmenting such pre-trained models to enforce additional physical constraints inference time, without requiring retraining.
3.2.1 Constraint Types in PDE Systems
Constraints in PDE-governed physical systems arise from a variety of sources, including bound- ary conditions, conservation laws 20 62], and physical admissibility requirements [63]. These constraints can be local (e-g., pointwise Dirichlet or Neumann conditions) or global (e.g., integral conservation of mass energy), may be linear or nonlinear in the solution field u(æ, t), and may not

---

Table 2: Summary of constraint types commonly encountered in PDE-based physical systems, categorized by their mathematical form and scope.
Constraint Type
Representative Form
Linearity
Dirichlet IC BC
Hu = Au Au bb=
Linear
Global mass conservation (periodic BCs)
Hu Jo t(2, C
Linear
Nonlinear conservation law
Hu d ,)
Nonlinear
Neumann or flux boundary condition
Hu (, - g(x,t)
Potentially nonlinear
Coupled or implicit constraints
Nonlinear spatial/temorra relationships
Nonlinear
be implicitly satisfied by the generated solution. Table summariize representative forms of com- monly encountered constraints, categorized by their mathematical structure. This taxonomy guides the design of appropriate enforcement mechanisms in generative modeling pipelines. A detailed description of the constraints, particularly those in our experiments, is provided in Appendix D
Algorithm 1 PCEM: Physics-Constrain Flow Matching
Require: Flow model vo (u, 7), constraint residual h(u), initial state ug, steps N, penalty X Ensure: Final state such that h(n)) 0
Ar f- 1/N, u f- u0
2: for k O....,N - d
3:
T - ÂT. T' + AT
4:
"1 ODESolve(u, 1o, 1,
5:
J Vh(un)
6:
Uproj + #n - J (JJ )-1((())))
7:
a, ODESollve Uproj. - (Uupooj uo): 1, 7')
8:
uy + arg ming ||u - ú, |l + (1 -' r')uo(u,e, +')
9:
u Ugy
10:
if [I(xw)] > then
+ arg ming ||u u1l1 s.t. h(u)= 0
12: return
3.3 PCFM: Physies-Constrained Flow Matching
Given a pretrained generative flow model ug(u, T), the flow dynamics define pushforward map 6, that transports samples uo To to solution-like outputs иı Φi (wo) Ti In our setting, constraints are imposed on the final sample W1, via a nonlinear function Hu h(w)) = 0, which we wish to enforce exactly during inference trajectory such that the terminal output 41 satisfics h(u1) = O, while remaining aligned with the learned flow ug. To achieve this, develop constraint-guice -uided sampling algorithm that interleaves lightweight constraint corrections with marginally consistent flow updates. The core procedure is outlined below.
3 Our goal is to generate samples u(T) along the generative
Forward Shooting and Projection. At each substep 7 + 8T, we first perform a forward solve (shooting) to the final flow time 7 = 1, 41 = ODESolve(u((T) uo, T, 1), typically using an inexpensive integrator (e.g., Euler with large step size) since high precision is not required for intermediate inference. We then apply single Gauss-Newton projection to softly align u to the constraint manifold.
r >=h((,) Vh(u)) Uproj u1 J JJ
which shifts in minimally to a linearized feasible point. As formalized in Proposition E. this projection step corresponds projection onto the tangent space of the constraint manifold #1, and recovers the exact solution when h is affine.
*While the usage of the term hard constraint typically implies that the constraint residual must be exactly zero, we use it synonymously with approximate hard constraint, referring to constraint satisfaction up to solver and machine numerical precision, common in related works (54 14
5

---

Reverse Update via OT Interpolant. To propagate this correction back to T' = T åT, we consider a reverse ODE solve u ODESolve (u(T), vgg 1. 1.7 ODEs are inherently reversible in theory, however in practice it may have (1) error the backward integration [64] mainly because of the flip of the signs of the eigenvalues of the Jacobian of Ua [47]. Gholami et al. showcases this
can be unstable for neural ODEs with large absolute value of the Lipschitz constant of ve 47 46].
However, robust options, such as implicit or adjoint methods [48] 65], are computationally prohibitive during inference. To alleviate these issues during inference, instead, we approximate the reverse flow using the Optimal Transport (OT)
displacement interpolant used during flow matching:
á, (uEoovee(upro) (oppro:-(uppp)) (upro) - 1,7)
(5)
This avoids instability and is justified by Proposition [3.1] which shows that the straight-line dis- placement approximates the true marginal flow the small-step limit. This perspective further supported by recent findings in Rectified Flows [36]. where generative paths become increasingly linear as flow resolution improves.
Proposition 3.1 (Reversibility under OT Displacement Interpolant). Let vo(u, T) be a pre-trained marginal velocity field learned via deterministic flow matching, with pushforward map or such that u1 = Φi (uo). for samples ug no and un TT. Suppose va valu, T) is Lipschitz continuous in both and r. Then under namerical integration with step size ÖT, the flow satisfies
vo(u(T),7) mn un ((87))
where p is the order of the integrator. Moreover, defining the OT displacement interpolant t(u) := "1 ug, the reverse update ü, ODE Uproj -ū(u), 1, produces numerically stable and solver-invarian approximation that is unconditionally reversible as Ô7 - 0.
Relaxed Constraint Correction. Due to discretization and the nonlinearity of h, the point ú,, may still incur residual error. We thus perform penalized correction:
Up arg min ||u - 1 + hiu + tor((,,)))) 7=1-7,
(6)
which encourages constraint satisfaction the extrapolated point while preserving local flow align- ment. As 8f - O, this step aligns with the OT-based marginal objective, and offers robustness in settings with coarse discretization, for e.g., necessitated by the need of lower number of function evaluations. We provide an ablation in Appendix L.3 highlighting its effect.
Final Projection and Runtime. We set Uy == uyy, and iterate. constraint residuals remain high at r = 1, we apply final full projection:
1 = arg min 1/1/ u, 12 subject to h(w') = 0.
'
To this end, we have developed custom, batched differentiable solver that performs these updates entirely in parallel: at each iteration assembles and inverts only an m X m Schur complement system (with dim h < n)), then applies Newton-style correction to the full state (see Appendix for the specific algorithm details). Despite the extra iterations, the overall cost remains ((n) per sample, since (i) each Schur solve is O(m") with << (ii) the backward flow solve in Eq. (5) is single
Euler step, and (iii) the relaxed- relaxed-penalty correction in Eq. (6) typically converges in 3-5 gradient steps. We also provide detailed runtime analysis in Appendix E..I, showing that the final projection step contributes only 1-3% of total time; most cost is due to sampling itself.
We note that ECI [12 emerges as special case of this framework when constraints are linear and non-overlapping. and the penalty is omitted (A 0) (see Appendix (F). Our method generalizes this by enabling strict enforcement of nonlinear and overlapping constraints while preserving generative consistency. A separate ablation study is also provided highlighting the necessity and role of custraint-guided sampling with consistent flow updates in Appendix
3.4 Reversibility and Interpolant Compatibility
Our method naturally extends to other interpolants, including variance preeercvir (VP) Score-Based Diffusion Models (SBDMs) 2 While OT interpolants yield straight-line paths, VP flows benefit from bounded Lipschitz constants due to Gaussian smoothing, making them more stable in reverse [6)] The key requirement access to a consistent displacement direction (e.g., #n un), which enables constraint enforcement via reversible approximations regardless of the underlying interpolant.

---

Table 3: Generative performance for zero-shot methods constrained PDEs with linear and nonlinear constraints, Heat and Navier-Stokes enforce global conservation laws (CL) as linear constraints, along with initial condition (IC) constraints. In contrast, Burgers and Reaction-Diffuse apply CL as nonlinear constraints along with IC or BC constraints. Lower values indicate better performance, and best results highlighted in bold.
Dataset
Metric
PCFM
ECI
DiffusionPDE
D-Flow
PDM
FFM
MMSE/ 10 10
0.241
0.697
4.49
1.97
0.45
4.56
SMSE/ 10
0.937
0.973
3.93
1.14
0.02
3.51
Heat Equation
CE (IC) 10
0
599
102
579
CE (CL) / 10
0
2.06
64.8
2.11
FPD
1.22
1.34
1.70
2.70
0.17
1.77
MMSE/ 10
4.59
5.23
17.4
12.21
16.5
SMSE/ 10
4.17
7.28
9.48
6.61
7.90
Navier-Stokes
CE (IC) 10
0
288
328
CE (CI.) 10-2
21.4
18.6
FPD
1.00
1.04
3.70
2.01
2.81
MMSE/ / 10 10
0.026
0.324
3.16
0.318
1.74
2.92
SMSE/ 10
0.583
0.060
2.54
6.86
0.43
2.54
Sraction-Diffusion IC
CE (IC) / 10-2
215
445
CE (CL) 10
0
6.00
3.82
29.7
3.87
FPD
15.7
136
28.3
24.9
MMSE/ 10 -2
0.335
0.359
5.42
0.224
11.8
4.86
SMSE 10 10-2
0.123
0.089
1.30
0.948
2.67
1.38
Burgers BC
CE (BC) 10
20.3
95.7
409
CE (CL) 100
15.7
6.20
15.0
6.91
FPD
0.292
0.307
25.9
1.44
6.41
24.7
MMSE/ / 10 10
0.052
10.0
14.3
9.97
153
13.7
SMSE 10-2
0.272
6.65
8.06
7.91
2.62
7.90
Burgers IC
CE (IC) 10
0
471
397
462
CE (CL)/10
205
6.22
8.66
6.91
FPD
0.101
1.31
35.8
22.1
99.7
33.5
Pretrained Poseidon requires spatially square inputs; bilinearly interpolated solution grids 128x 128, which may introduce artifacts evaluation. D-Flow and PDM results are omitted due to numerical instabilities.
Pretrained Poseidon requires spatially square inputs; bilinearly interpolated solution grids 128x 128, which may introduce artifacts FPD evaluation. D-Flow results are omitted due to numerical instabilities.
3.4.1 Tradeoff Between Flow Steps and Relaxed Constraint Correction
In PCFM, the number of flow steps and the strength of the relaxed constraint correction jointly determine the quality of the final solution. Using fewer steps speeds up inference but increases numerical error, potentially violating constraints at intermediate states. In such regimes, the relaxed correction term in Equation (6) provides mechanism compensate for deviation from the constraint manifold. Conversely, with sufficiently many flow steps, the penalty term becomes less critical, and the dynamics naturally stay closer the OT path (36). We present detailed ablation results in Appendix L.3 showing how PCFM adapts across this tradeoff.
Experiments and Results
In this section, we highlight the key results and compare our approach with representative baselines. For every problem, we construct PDE numerical solution dataset with two degrees of freedom by varying initial (IC) and boundary conditions (BC) (see Appendix B). and pre-train an uncondi- tional FFM model [II] on this dataset (see Appendix C). To evaluate generative performance and constraint satisfaction, we focus on generating solution subsets constrained on selected held-out IC or BC, aiming to guide the pretrained model toward the corresponding solution subset. Fur- thermore, we incorporate physical constraints, specifically global mass conservation, via PCFM and adapt other baseline methods, where possible, through their sampling frameworks evaluate
7

---

their performance on constraint satisfaction tasks. We open-source our code: a Python implemen- tation at https: wwwpppfpenggei//////ooooooocpc //CMMM and an experimental Julia implementation at https://github. comm/utkarrshs3/PPPN... j1.
For variety of PDEs with different constraint requirements, ranging from easier linear constraints to harder nonlinear constraints, we evaluate different sampling methods using the same pretrained FFM model. The specific linear and nonlinear constraints used in our PCFM approach are outlined in Appendix D while the parameters and set-up for other methods are outlined in Appendix H. To evaluate the similarity between the ground truth and generated solution distributions, we follow Kerrigan et al. [TT] and Cheng et al. [1] to compute the pointwise mean squared errors of the mean (MMSE) and standard deviation (SMSE), along with the Fréchet Poseidon Distance (FPD). The FPD quantifies distributional similarity in feature space using a pretrained PDE foundation model (Poseidon) [67]. analogous to Fréchet Inception Distance (FID) used in image generation [68] Importantly, we evaluate the constraint errors by taking the l norm of the residuals for initial condition (IC), boundary condition (BC), and global mass conservation (CL.) over time f, averaged across all N generated samples: CE(+) = N Σ == Rx (1
where * E {IC, BC, CL}
4.1 Linear Hard Constraints
We first focus on linear hard constraints by considering the 2-D Heat equation (u(2r, ()) and the 3-D Navier-Stokes (u(x,, y, ()) equation with periodic boundary conditions. The global conservation laws simplify to a linear integral over the solution field, where the conserved mass quantity reduces to sum or mean of the solution across the spatial domain (see Appendix D). We aim to constrain the generation to satisfy a selected IC and global mass conservation. In Table PCFM outperforms all the compared methods in the MMSE, SMSE, and FPD, achieving machine-leve precision for constraints on both IC and global mass conservation. Due to the simplicity of the linear constraints, ECI is also capable of achieving constraint satisfaction in generated samples.
PCFM further enables us to improve the performance by applying relaxed constraint corrections (Equation 6). achieving higher solution fidelity with fewer flow matching steps. For high-dimensional settings such as Navier-Stokes, we adopt a stochastic interpolant perspective by randomizing "o across batches [37]. This enhances alignment with the OT displacement path while preserving approximate hard constraint satisfaction at generation time.
4.2 Nonlinear Hard Constraints
4.2.1 Nonlinear Global and Boundary Constraints
Reaction-Diffusion (nonlinear mass & Neumann flux). We first consider a reaction-diffusion PDE with nonlinear mass conservation and Neumann boundary fluxes. When constrained on the fixed IC and nonlinear mass, PCFM achieves the lowest constraint errors and best generation fidelity (Table Figure 2 ECI meets the IC constraint by exact value enforcement in their correction step but its framework cannot enforce nonlinear mass conservation. DiffusionPDDF and D-Flow, despite using a combined loss function on IC and PINN-loss, likewise fail to enforce both hard constraints simultaneously. In contrast, PCFM satisfies IC and mass conservation to machine precision for all t (Figure 2). which in turn enforces Neumann boundary fluxes and improves overall solution quality.
Burger's s Equation (nonlinear mass & Dirichlet BCs). Constraining Dirichlet BC at x 0 and zero-flux Neumann BC at x = 1, generated samples should satisfy both BC and nonlinear mass conservation constraints for all t. PCFM achieves constraint satisfaction while matching similar fidelity of other methods (Figure 5).
4.2.2 Nonlinear Local Dynamics and Shock Constraints
We demonstrate PCFM's ability tackle more challenging task of enforcing global nonlinear and local dynamical constraints, achieving improved generation quality while enforcing physical consistencies. Specifically, for Burgers' equation constrained on ICs, PCFM outperforms all baselines across metrics while satisfying both IC and mass conservation. In addition to global constraints, we incorporate unrolled local flux collocation points in the residual (see Appendix for further details). Projecting intermediate in to the constraint manifold helps capture the localized shock structure (see top left corner of Figure improving solution fidelity to match the shock dynamics in Burgers. In contrast, other methods cannot satisfy hard constraints nor capture the shock dynamics.
8

---

Ground Truth
PCFM
FFM
DFlow
EC
DiffusionFIE
PDM
Generated
solution
mean
PCFM
Dif.sonPDE
Mass
residual
mean
Figure 2: Comparison of mean + std. of mass residuals across samples. generated solutions and mass conservation errors for the Reaction-Diffusinn problem with IC fixed. By enforcing both IC and nonlinear mass conservation constraints, PCFM improves quality of generated solutions while satisfying both constraints exactly.
Generated
solution mean
PCPM
FFM
DMfuserDD
POM
Mass residual
Figure 3: Comparison of mean generated solutions and mass conservation errors for the Burger's problem with IC fixed. By enforcing nonlinear conservation constraints via PCFM, our method captures the Burgers' shock phenomenon, ensures global mass conservation in the generated solution, while improving solution quality. Shaded bands show + std. of mass residuals across samples.
4.3 Enhancing Fidelity through Additional Physical Constraints
In our ablation study (Figure investigate the effect of imposing more constraints, specifically the effect of constrained collocation points, in addition to the IC and mass conservation constraints in the Burgers' problem. We show that adding more constrained collocation points improves generation quality without worsening the satisfaction of IC and mass constraints (see Appendix D for constraint details). Interestingly, while stacking multiple constraints increases computational Gauss-Newton projection, we find that complementary constraints (local flux, IC, and global mass) can improve performance, unlike soft-constraint approaches such PINNs, where adding more constraints can degrade performance due to competing objectives 26 25
in each
However, in cases of conflicting constraints, tradeoffs will arise that require careful balance of their effects on generation fidelity and physical consistency. Indeed, the ability to chain different hard constraints on a pretrained generative model makes PCFM flexible and practical for diverse applications.
We also explored total variation diminishing (TVD) constraints to promote smoother solution profiles, demonstrating on the heat equation that TVD improves both smoothness and fidelity while preserving IC and global mass constraints (see Appendix D). This highlights PCFM's flexibility to incorporate multiple variety of constraints enhance generative quality.
Conclusion
We presented PCFK, a zero-shot inference framework for enforcing arbitrary nonlinear equality constraints in pretrained flow-based generative models. PCFM combines flow-based integration, tangent-space projection, and relaxed penalty corrections to strictly enforce constraints while stay- ing consistent with the learned generative trajectory. Our method supports both local and global constraints, including conservation laws and consistency constraints, without requiring retraining or architectural changes. Empirically, PCFM outperforms state-of-the-at baselines across diverse PDEs,

---

MMSE
SMSE
IC Residual
Mass Residual
30
10ah
10
10
10
100 200
300
400
500
100
200
300
4.0
500
100
200
300
400
500
100
200
300
400
300
Constrained collocation points
Constralned collocation points
Constralined collocation points
Constrained collocation points
Flow matching steps
20 -
100
200
Figure 4: Increasing the number of constraints (constraint collocation points) can improve solution fidelity while maintaining strong satisfaction of other constraints (IC and global mass conservation), demonstrating the ability of PCFM to handle chaining of multiple constraints.
including systems with shocks and nonlinear dynamics, achieving lower error and exact constraint satisfaction. Notably, we find that enforcing additional complementary constraints improves gen- eration quality, contrary to common limitations observed in soft-penalty methods such as PINNs 25 [26]. While PCFM currently focuses equality constraints, extending it to inequality constraints and leveraging structure in constraint Jacobians are promising directions for improving scalability 669 [5]] Our work offers principled approach for strictly enforcing physical feasibility in generative models, with broad impact on scientific simulation and design.
Acknowledgments
We acknowledge the MIT SuperCloud and Lincoln Laboratory Supercomputing Center for providing HPC resources. We also acknowledge the Delta GPU system at the National Center for Supercomput- ing Applications (NCSA), supported by the National Science Foundation and the ACCESS program. P.C. is supported by the Regenerative Energy-Efficient Manufacturing of Thermoset Polymeric Materials (REMAT) Energy Frontier Research Center, funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under award DE-SC0023457.
This material is based upon work supported by the U.S. National Science Foundation under award Nos CNS-2346520, PHY-2028 15, RISE-2425761, DMS-23251 184,, OAC-2103800 and OSI-2029670, by the Defense Advanced Research Projects Agency (DARPA) under Agreement No. HR001 12490488, by the Department of Energy, National Nuclear Security Administration under Award Number DE- NA00039655 and by the United States Air Force Research Laboratory under Cooperative Agreement Neither the United States Government nor any agency thereof, nor
Number 148750-19-2-10000 any of their employees, makes any warranty, express implied, assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute imply endorsement, recommendation, favoring by the United States Government any agency thereof. The views and opinions of authors expressed herein not necessarily state or reflect those the United States Government or any agency thereof." The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the United States Air Force or the U.S. Government.
References
[1] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances neural information processing systems, 33:6840-668 6851, 2020.
[2] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011. 2020.
[3] Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. Advances in neural information processing systems, 34:8780-8794, 2021.
10

---

[4] Yaron Lipman, Ricky TQ Chen, Heli Ben- Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:2210.02747, 2022.
[5] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al. Language models are few-shot learners. Advances in neural information processing systems, 33:1877--190 1901, 2020.
[6] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
[7] John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ron- neberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Zídek, Anna Potapenko, al. Highly accurate protein structure prediction with alphafold. nature, 596(7873):333-58,, 2021.
[8] Gabriele Corso, Hannes Stärk, Bowen Jing, Regina Barzilay, and Tommi Jaakkola. Diffdock: Diffusion steps, twists, and turns for molecular docking. arXiv preprint 22 2022.
[9] Juno Nam, Sulin Liu, Gavin Winter, KyuJung Jun, Soojung Yang, and Rafael Komeez-Bombarelli Flow matching for accelerated simulation of atomic transport in materials, 2025. URL https
//arxiv www.:orggbbb//eee
01464
[10] Jiahe Huang, Guandao Yang, Zichen Wang, and Jeong Joon Park. Diffusionpde: Generative pde-solving under partial observation. arXiv preprint urXiv:2406.17763. 2024.
[11] Gavin Kerrigan, Giosue Migliorini, and Padhraic Smyth. Functional flow matching. arXiv preprint arXiv:2305. 17209, 2023.
[12] C. Cheng, B. Han, D.C. Maddix, A.F. Ansari, A. Stuart, M.W. Mahoney, and Y. Wang. Gradient- free generation for hard-constrained systems. In International Conference Learning Repre- sentations, 2025.
[13] Gefan Yang and Stefan Sommer. A denoising diffusion model for fluid field prediction. arXiv preprint 1616,, 2023.
[14] Derek Hansen, Danielle C. Maddix, Shima Alizadeh, Gaurav Gupta, and Michael W. Mahoney. Learning physical models that can respect conservation laws. In Proceedings of the 40th International Conference Machine Learning, volume 202, pages 22469-12518 PMLR, 2023.
[15] Taco Cohen and Max Welling. Group equivariant convolutional networks. In International conference machine learning, pages 2999. PMLR, 2016.
[16] Ricky TQ Chen and Yaron Lipman. Flow matching general geometries. arXiv preprint arXiv:2302.03660, 2023.
[17] Mevlana Gemici, Danilo Rezende, and Shakir Mohamed. Normalizing flows on riemannian manifolds. arXiv preprint arXiv:1611.02300 2016.
[18] Emile Mathieu and Maximilian Nickel. Riemannian continuous normalizing flows. Advances in Neural Information Processing Systems, 33:2503-2515, 2020.
[19] Lawrence C Evans. Partial differential equations, volume 19. American Mathematical Society,
2022.
[20] Randall J. LeVeque. Numerical Methods for Conservation Laws. Lectures in mathematics ETH Zürich. Birkhäuser Verlag, 1990.
[21] Emmy Noether. Invariant variation problems. Transport theory and statistical physics, 1(3): 186-207, 1971.
[22] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
11

---

[23] 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, 2021.
[24] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, (()):218-22)) 2021.
[25] Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies physics-informed neural networks. SIAM Journal on Scientific Computing. 43 (5):A3055-A3081, 2021.
[26] Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in thysics-informed neural networks. In Advances in neural information processing systems, volume 34, pages 26548-26566 2021.
[27] Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449:1 2022.
[28] Jack Richter-Powell, Yaron Lipman, and Ricky T. Q. Chen. Neural conservation laws: A divergence-free perspective. arXiv preprint arXiv:2210.0772 2022.
[29] Takashi Matsubara and Takaharu Yaguchi. Finde: Neural differential equations for finding and preserving invariant quantities. In International Conference on Learning Representations, 2023.
[30] S. Chandra Mouli, Danielle C. Maddix, Shima Alizadeh, Gaurav Gupta, Andrew Stuart, Michael W. Mahoney, and Yuyang Wang. Using Uncertainty Quantification to Character- ize and Improve Out-of-Domain Learning for PDEs. In Proceedings of the 40th International Conference Machine Learning, volume 235, pages 36372-36418. PMLR, 2024.
[31] Geoffrey Négiar, Michael W. Mahoney, and Aditi S. Krishnapriyan. Learning differentiable solvers for systems with hard constraints. International Conference on Learning Representa- tions, 2023.
[32] Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Advances in neural information processing systems, volume 32, 2019.
[33] Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXiv preprint arXiv:2003.04630, 2020.
[34] Lyle Regenwetter, Giorgio Giannone, Akash Srivastava, Dan Gutfreund, and Faez Ahmed. Constraining generative models for engineering design with negative data. Transactions on Machine Learning Research, 2024.
[35] Heli Ben-Hamu, Omri Puny, Itai Gat, Brian Karrer, Uriel Singer, and Yaron Lipman. D-flow: Differentiating through flows for controlled generation. arXiv preprint arXiv:2402. 14017, 2024.
[36] Xingchao Liu, Chengyue Gong, and Qiang Liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. arXiv preprint arXiv:2209.0300 2022.
[37] Michael S Albergo, Nicholas M Boffi, and Eric Vanden- Fijnden. Stochastic interpolants: A unifying framework for flows and diffusions. arXiv preprint arXivv2303.08897, 2023.
[38] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[39] Michael S Albergo and Eric Vanden-Eijnden. Building normalizing flows with stochastic interpolants. arXiv preprint arXiv:2209.15571, 2022.
[40] Jae Hyun Lim, Nikola B Kovachki, Ricardo Baptista, Christopher Beckham, Kamyar Az- izzadenesheli, Jean Kossaifi, Vikram Voleti, Jiaming Song, Karsten Kreis, Jan Kautz, al. Score-based diffusion models in function space. arXiv preprint arXiv:2302.07400, 2023.
[41] Dule Shu, Zijie Li, and Amir Barati Farimani. A physics-informed diffusion model for high- 2023.
fidelity flow field reconstruction. Journal of Computational Physics,
12

---

[42] Bowen Song, Soo Min Kwon, Zecheng Zhang, Xinyu Hu, Qing Qu, and Liyue Shen. Solv- ing inverse problems with latent diffusion models via hard data consistency. arXiv preprint arXiv:2307.08123, 2023.
[43] Xingchao Liu, Lemeng Wu, Shujian Zhang, Chengyue Gong, Wei Ping, and Qiang Liu. Flow- grad: Controlling the output of generative odes with gradients. In Proceedings of the IEEE/CVF Conference Computer Vision and Pattern Recognition, pages 14335-24344 2023.
[44] Hengkang Wang, Xu Zhang, Taihui Li, Yuxiang Wan, Tiancong Chen, and Ju Sun. Dm- plug: A plug-in method for solving inverse problems with diffusion models. arXiv preprint arXiv:2405. 16744, 2024.
[45] Hyungjin Chung, Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye. Improving diffusion models for inverse problems using manifold constraints. Advances in Neural Information Processing Systems, 35:25683 222296666,,,,,,, 2022.
[46] Amir Gholami, Kurt Keutzer, and George Biros. Anode: Unconditionally accurate memory- efficient gradients for neural odes. arXiv preprint 17Xiv:1902.10298, 2019.
[47] Suyong Kim, Weiqi Ji, Sili Deng, Yingbo Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Intendisciplinary Journal of Nonlinear Science, 31)),, 2021.
[48] Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint arXiv:2001.04388,, 2020.
[49] Jacob K Christopher. Stephen Baek, and Nando Fioretto. Constrained synthesis with projected diffusion models. Advances in Neural Information Processing Systems, 37:89307-89333, 2024.
[50] Dimitri Bertsekas. Nonlinear programming. Journal of the Operational Research Society, 48 (3:33443334 1997.
[51] George Em Kamiaadkii, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
[52] Zhou Lu, Hongming Pu, Feicheng Wang, Zhiqiang Hu, and Liwei Wang. The expressive power of neural networks: view from the width. In Advances in Neural Information Processing Systems, volume 30, 2017.
[53] Pengfei Cai, Sulin Liu, Qibang Liu, Philippe Geubelle, and Rafael Gomez-Bombarelli. Towards long rollout of neural operators with local attention and flow sattchin-insspiee correction: An example in frontal polymerization pdes. NeurIPS 2024 Workshop on Machine Learning and the Physical Sciences, 2024.
[54] Priya L Donti, David Rolnick, and J Zico Kolter. Dc3: A learning method for optimization with hard constraints. In International Conference on Learning Reprexentations, 2021.
[55] Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. In Advances in neural information processing systems, volume 32, 2019.
[56] Utkarsh Utkarsh, Danielle C Maddix, Ruijun Ma, Michael Mahoney, and Yuyang Wang. End-to-end probabilistic framework for learning with hard constraints. arXiv preprint arXiv:2506.07000 2025.
[57] Aaron Lou, Derek Lim, Isay Katsman, Leo Huang, Qingxuan Jiang, Ser Nam Lim, and Christopher M De Sa. Neural manifold ordinary differential equations. Advances in Neural Information Processing Systems, 33: 13:17758817558 2020.
[58] Muhammad Firmansyah Kasim and Yi Heng Lim. Constants of motion network. In Advances in Neural Information Processing Systems, volume 35, pages -93305, 2022.
[59] Uri M Ascher and Linda R Petzold. Computer methods for ordinary differential equations and infferettialalgebbrt equations. SIAM, 1998.
13

---

[60] Linda Petzold. Differential/algebriic equations are not ode's. s. SIAM Journal on Scientific and Statistical Computing, 3(3):367-384, 1982.
[61] Ernst Hairer, Marlis Hochbruck, Arieh Iserles, and Christian Lubich. Geometric numerical integration. Oberwolfach Reports, 3(1):805-882, 2006.
[62] Randall J Le Finite volume methods for hyperbolic problems, volume 31. Cambridge university press, 2002.
[63] Leonid Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(1-4):259-26) 1992.
[64] Gerhard Wanner and Ernst Hairer. Solving ordinary differential equations II, volume 375. Springer Berlin Heidelberg New York, 1996.
[65] Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.02433. 2022.
[66] Joe Benton, George Deligiannidis, and Arnaud Doucet. Error bounds for flow matching methods. arXiv preprint arXiv:2305. /6860, 2023.
[67] Maximilian Herde, Bogdan Raonić, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Em- manuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for pdes, 2024.
arxiv. arg/abs/24055 19101
[68] Jae Hyun Lim, Nikola B. Kovachki, Ricardo Baptista, Christopher Beckham, Kamyar Azizzade- nesheli, Jean Kossaifi, Vikram Voleti, Jiaming Song. Karsten Kreis, Jan Kautz, al. Score-based diffusion models in function space, 2023. httpsstt::///xs. org/ hrrg/ab//302.... 2302 07400.
[69] Stephen P Boyd and Lieven Vandenberghe. Convex optimization, Cambridge university press, 2004.
[70] Carl T Kelley. Iterative methods for linear and nonlinear equations. SIAM, 1995.
[71] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, and Luca Antiga. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
[72] Utkarsh Utkarsh, Valentin Churavy, Yingbo Ma, Tim Besard, Prakitr Srisuma, Tim Gymnich, Adam R Gerlach, Alan Edelman, George Barbastathis, Richard D) Braatz, et al. Automated translation and accelerated solving of differential equations on multiple gpu platforms. Computer Methods in Applied Mechanics and Engineering, 419:1 2024.
[73] Kai Fukami, Koji Fukagata, and Kunihiko Taira. Super-resolution reconstruction of turbulent flows with machine learning. Journal of Fluid Mechanics, 870: 120, 2019.
[74] Stephan Rasp, Peter D Dueben, Sebastian Scher, Jonathan A Weyn, Soukayna Mouatadid, and Nils Thuerey. Weatherbench: a benchmark data set for data-driven weather forecasting. Journal of Advances in Modeling Earth Systems, 12(1 1)e2020MS000220, 2020.
14

---

A Proof of Constraint Satisfaction
We now show that, under standard regularity assumptions. PCFM (see Algorithm sample リ satisfying the hard constraint huu1) = 0 to numerical precision.
D returns a final
Theorem A.I (Exact Constraint Enforcement). Let h RRn _ TR. be twice continuously differen- tiable constraint function, and suppose the Jacobian (u) == Vh(u) Rmxx has full row rank
m S. n in a neighborhood of the constraint manifold M == {u R2 h(u) = 0}. Then the final sample in produced by PCFM satisfies
h((1)) machine precision, provided the final projection step in Algorithm is solved using sufficiently many Newton-Schur iterations.
Proof. We consider the final correction step in Algorithm 1:
u1 == arg min ||u' - u1ll s.t. h(u')= =
which seeks to project the current sample un (obtained after flow integration and relaxed penalty correction) onto the feasible manifold.
Let (0) ui denote the initial guess to the Newton solver. The constrained root-finding problem h(u) = 0 is nonlinear, but under the full-rank Jacobian assumption, the implicit function theorem ensures locally unique root u² E EM near (0), We solve the nonlinear system using the following Newton-Schur iteration:
4 (++1) = (T) ( (a( (h(ul()) + Jn (n - =(())) (A.1)
This corresponds to the batched Schur-complenen Newton update implemented in differentiable solver. Standard convergence theory (e.g.. Theorem 10.1 in (70) guarantees that, for sufficiently small [[b(u())]] the iterates (A) converge quadratically to u². Hence, for any tolerance e > 0, there exists finite K such that
[lh(u()))] <E
Since we continue iterations until [[ < Ewl (with Eul typically set to machine precision), we conclude
h(m)) (k) ~ O.
Finally, we verify that the initial guess (0) un is within the basin of convergence. Each prior step of PCFM includes: (i) Gauss-Newtte projection, (ii) an OT-based backward solve, and (iii) relaxed penalty correction. These ensure the incoming residual [b(() |I is already small (typically < 10 30- Combined with the well-conditioned Jacobian Jn. we ensure convergence of the final Newton-Schur projection.
B PDE Dataset Details
We evaluate our method on four representative PDE systems commonly studied scientific machine learning: the Heat Equation, Navier-Stokes Equation, Reaction-Diffusion Equation, and Burgers' Equation. These problems cover linear and nonlinear dynamics, have smooth or discontinuous solutions, and include a variety of initial and boundary condition types.
B.1 Heat Equation
We consider the one-dimensional heat equation with periodic boundary conditions, following the setting in Hansen et al. [14]
u, atl
E [0, t [0, 1],
(8)
with the sinusoidal initial condition
u(z, 0) sin(x sp)
(9)
15

---

and periodic boundary conditions w(0, t) = u(2Ttt), t). The analytical solution is given by u( t) exp(-at) sin(x: + (). To pretrain our FFM model, we construct the dataset by varying the diffusion coefficient a 4(1,5) and phase yp 4(0,, r). All solutions are on 100 x 100 spatiotemporal grid. During constrained sampling, we restrict to the same initial condition by fixing 4 = 7/4. The global mass conservation constraint is enforced via
2x
u(x, t) =0,
Vt E [0. 1].
(10)
B.2 Navier-Stokes Equation
We consider the 2D incompressible Navier-Stokes (NS) equation in vorticity form with periodic boundary conditions, following Li et al. [23]:
0,w(r, + u(2, t)
Vw((x t) 1/,w(z, f(x),
(11)
V uu(( () = 0,
(12)
where w = V x u denotes the vorticity and y = 10-3 3 is the viscosity. The initial vorticity Wo sampled from Gaussian random field, and the forcing function is defined as: f(z) 0.1v2 sin (2 (a1 + x2) + +), o ~ (00 /2). We solve the NS system using Crank-Nicolson spectral solver on a 64 x 64 periodic grid, recording 50 uniformly spaced temporal snapshots over the interval T = 49. For training, we generate 10000 simulations by sampling 100 random initial vorticities and 100 forcing phases. For test dataset, sample an additional 10 vorticities and 100 forces, yielding 1000 solutions. We fix to initial vorticity for comparison during sampling.
Here, we prove the global mass conservation for this NS problem set-up. Under periodic boundary conditions, global vorticity is conserved. Specifically, integrating the vorticity equation over the spatial domain S2 and applying the divergence theorem yields:
J. w(z,t) / vAw((z,t) + f(x) u (w(x, dr 0,
since the divergence and Laplacian terms vanish due to periodicity and the incompressibility condition V M 0. Thus, the total vorticity Jo wr(ar, , 1) de remains constant for all
B.3 Reaction-Diffusiae Equation
We consider a nonlinear ID praction-diffusiio equation with Neumann boundary conditions:
We pp((1 u) v0=u, u,
(13)
on the domain x E [0, 1J, e [o, I], where (p. and right ends are specified IL and ØR. respectively. The initial condition u(x; is sampled from a randomized combination of sinusoidal and localized bump functions (see Appendix D for full specification). The training dataset is constructed by pairing 80 initial conditions with 80 boundary conditions, producing 6400 PDE solutions, all discretized on (nz, = (128, 100) space-time grid. We use semi-implicit finite difference solver with CFL-controlled time stepping. During constrained sampling, we fix one initial condition and vary the boundary fluxes to enforce constraint satisfaction.
0.005). The boundary fluxes at the left
As this system involves nonlinear source and flux terms, global conservation is no longer trivial, Integrating both sides of the PDE over the domain yields:
dt / t(zxt)dd do = .f u(1 u) dx + gL( gR((),
d d
(14)
which depends nonlinearly on the state u(z, t) and boundary fluxes gL. 9R. be tracked through both the nonlinear reaction term and boundary-driven transport, making hard constraint enforcement via our framework necessary.
gR. Thus, conservation must
B.4 Inviscid Burgers' Equation
We consider the 1D inviscid Burgers' equation with mixed Dirichlet and Neumann boundary condi- tions:
1 U + 2 (12)) = x (0, 1],
t [0, 1).
(15)
16

---

with initial condition defined as a smoothed step function centered at random location Ploc ~ (.002,0.8)
u(2 0) = 1 exp
= 0.02.
(16)
We apply a Dirichlet condition on the left, u(0, the with the ((, 1), and Neumann condition (zero-gradient) on the right. Solutions are computed using a Godunov finite volume method (nzz, nl) (101, 101) grid. For pretraining the FFM model, we vary both initial and boundary conditions, generating 80 variants each to create 6400 solutions in the training set. For constrained sampling. use two configurations: (I) fixing the initial condition and varying the boundary condition, and (2) fixing the boundary condition and varying the initial condition.
The conservation constraint is nonlinear and particularly sensitive due shock formation. Integrating the PDE yields
u(z, t) dr = [1, t)² (0. / u(x,
dt
so global conservation requires u(1,t)2 = u(0, t)?, making the constraint highly sensitive to coundayy-induced discontinuities.
C FFM Model Pretraining
For all PDE benchmarks, we employ the functional flow matching (FFM) 11
framework and
parameterize the underlying time-dependent vector field with Fourier neural operator (FNO) backbone (23). The FNO encoder takes the input the concatenation of the current state th,, sinusoidal positional embedding for the spatial grid and Fourier time embedding. For each benchmark, we adopt the following training scheme unless otherwise specified:
Optimizer: Adam optimizer with learning rate 3 x 10 81 = 0.9, 32 = 0.999, and no weight decay.
Learning rate scheduler: cattuce-on-platean scheduler with a factor of 0.5, patience of 10 validation steps, and minimum learning rate of 1 10-4
Batch size: 256 for ID problems (Heat, Reaction-Diffusionn Burgers) and 24 for the 2D Navier-Stokes problem.
For Heat, follow Cheng al. [12] to train 5000 analytical solutions of the heat equation. For and Burgers, we train each model on 6400 numerical PDE solutions. For these problems, the FNO is configured with 4 Fourier layers with 32 Fourier modes for both spatial and temporal dimensions, 64 hidden channels, and 256 projection channels. For Navier-Stokes equation, train 3D FFM 10, 000 numerical solutions. The FNO has Fourier layers with 16 Fourier modes per dimension, 32 hidden channels, and 256 projection channels, following Cheng et al. [12] FFM models for Heat, Reaction- NVIDIA V100 GPU and the Navier-Stokes model is trained for 500, 000 steps on 4 NVIDIA A100 GPUs.
and Burgers are each trained over 20, 000 steps on
Following Kerrigan et al. []], all FFM models use a Gaussian prior noise P(ug). For ID problems, adopt Matern Gaussian Process prior with smoothness parameter 0.5, kernel length 0.001, and variance 1.0, implemented using GPyTorch. For Navier-Stokes (2D), we adopt standard Gaussian white noise the prior, P(up) (/00 T). Importantly, PCFM is applied post-hoc at inference time, allowing us to reuse the pretrained unconditional FFM models across constraint configurations.
D Constraints Enforcement via PCFM
We focus physical and geometric constraints of the form Hu 0, which must be satisfied as equality conditions on the solution field u. First, we describe common equality constraints encountered in solving PDEs. Next, describe specific constraint setup for our PDE datasets.
17

---

D.1 Linear Equality Constraints
Dirichlet Initial and Boundary Conditions g(2,t) 0Q x [0, define fixed-value conditions that are linear in u. These can be encoded in the general form
Constraints such as u(x, 0) = au(()) or u( t) =
Hu Au b= 0,
(17)
where A is a sampling or interpolation operator applied to and is the prescribed target value.
Linear Global Conservation Laws (Periodic Systems) In systems with periodic boundary condi- tions, additive invariants such as total mass are often linearly conserved. A representative example
Hu [ u(xr 1) de C 0,
(18)
where E R is the conserved quantity. Such constraints can be exactly enforced via analytically
derived projections [14), geometric integrators [6]]. and numerical optimization (69 50
D.2 Nonlinear Equality Constraints
Many physical systems exhibit global conservation laws where the conserved quantity depends nonlinearly on u. A representative form is:
d Hu dt / p(u(x,t)) t)) O.
(19)
where p(u) denotes a nonlinear density, such as total energy or entropy. While such constraints often not conserved exactly by standard numerical integrators, such as symplectic or variational integrators [6]. can conserve invariants approximately over long time horizons. In setting, however, we seek enforce such nonlinear conservation laws exactly at inference time and moreover only at the final denoised sample. Thus, we rely on applying iterative projection methods such as Newton's method for constrained least squares (69 to satisfy Equation (19) up to numerical tolerance.
methods,
Neumann-Type Type Boundary Constraints
Boundary constraints involving fluxes, such as
Hu 0,,((, t) gxx, = 0,
(20)
are typically linear under standard coni-discretizations, as the derivative operator is linear. However, these constraints can become effectively nonlinear in practice when implemented with upwinding or other nonlinear dlx-diisrettization schemes, particularly in hyperbolic PDEs (62]. Additionally, if the prescribed flux g depends nonlinearly on u. the constraint becomes explicitly nonlinear. both
cases, we handle them using differentiable correction procedures projection-basee -ased updates within our framework.
In our PCFM framework, the residual operator H(u) specifies the equality constraints to be enforced on the generated solution u. We enforce these constraints the generated solution state in our PCFM projection operator, ensuring that the final solution strictly satisfies the required physical constraints without retraining the underlying flow model. We summarize the specific constraint formulations used for each PDE task below.
D.3 Heat Equation (IC and Mass Conservation)
The residual is concatenation of two linear constraints:
[
u(z, t 0) uic (x) J u(, L) do ʃ u(x, / = 0) dx
H (u)
(21)
where the first term enforces the initial condition and the second term ensures mass conservation over time.
18

---

D.4 Navier-Stokes Equation (IC and Mass Conservation)
For the 2D Havier-Stokee problem, we impose the initial vorticity and the global mass conservation constraint over the spatial domain:
u(I, v, u(2, dx dy -
0) unc (2c y) u(2, y, 0) dæ dy
HNs(u) =
/ʃ
(22)
where uzz y, t) represents the vorticity field over spatial coordinates (æ, y) and time t. The mass conservation term enforces that the total vorticity integrated over the spatial domain remains consistent with the initial condition over all time steps.
D.5 Dreactio-DDiffssion (IC and Nonlinear Mass Conservation)
We enforce both the initial condition and nonlinear mass conservation which accounts for reaction and boundary flux terms:
u(x, 0) uuc(x) HRD(U)= m(t)) - (m(0) + fo p(u) dt + [0(9 gR) dt
(23)
where m(t) f u(xx t) dr is the total mass, p(u) = pr( - u) is the nonlinear reaction source term, and JL, IR are the Neumann boundary fluxes.
D.6 Burgers (BC and Mass Conservation)
We enforce both the boundary conditions (Dirichlet and Neumann) and nonlinear mass conservation: 0,t) -
u(r HBurgens = u(x = -1,t) u(x -2 t) u(x t) dr u(x,t= = dæ
(24)
[
where ur, the Dirichlet boundary value at the left boundary, and the Neumann BC enforces zero-gradient at the right boundary.
D.7 Burgers (IC, Mass Conservation, and Local Flux)
We impose three complementary constraints for the Burgers equation: (i) the initial condition, (ii) global nonlinear mass conservation, and (iii) a sequence of local conservation updates based on Godunov's 's flux method. Specifically, the residual is formulated
u(x, t- = 0) Mc(()
R (k)
H HBurgas--Ic(14
R
(25)
f
Ju(a, dr
u(z, t 0) dx
where (k) (u) applies & unrolled discrete updates (collocation points) based on Godunov flux:
min ( fuh)
if if < uR if uf > ug and ur, > uR and un < 0
u(u
1 !
0
(26)
We apply & E {L....,5) unrolled steps incrementally enforce local conservation dynamics alongside global mass conservation and initial condition satisfaction.
E Details of the PCFM
E.1 Runtime Analysis and Overhead
Assessing runtime and memory overhead essential for evaluating the practical feasibility of PCFM. We therefore report detailed runtime and memory comparisons in Table and Table 5 across all
19

---

baselines using consistent batch sizes (8 for Kavier-Stokee and 32 for other datasets) and 100 flow- matching steps on single NVIDIA V100 GPU. To ensure fair and interpretable comparison. all experiments are conducted using À 0 in Equation 6 The choice of À can significantly affect convergence behavior and computational cost, as larger penalty weights introduce additional correction steps (a separate ablation is available in Appendix L.33). However, we observe that X = 0 sufficient to achieve state-of-the- performance in most cases, while also isolating the intrinsic cost of our projection method with greater clarity.
PCFM introduces modest overhead relative unconstrained baselines such as vanilla FFM, ow- ing to the additional projection step for enforcing constraints. Nevertheless, it achieves a 4-6x speedup over D-Flow [35, which requires costly backpropagation through the entire ODE unroll. DiffusionPDE [TO] runs faster but fails to enforce nonlinear constraints robustly, particularly for problems such as Burgers and Reaction-Diffusion ECI exhíbits runtime growth with the number of mixing steps but similarly struggles with nonlinear coupled constraints. For linear systems such as Heat and Navier-Stokes, the projection converges rapidly, resulting in comparable runtimes to PDM 49 Overall, PCFM incurs small yet necessary computational overhead to achieve robust constraint enforcement, maintaining competitive runtime efficiency while providing strict constraint satisfaction.
Table 4: Runtime and memory comparison across methods. Reported values correspond per-sample inference on a single NVIDIA V100 32 GB GPU.
Dataset
Metric
PCFM
ECI
DiffusionPDE
D-Flow
PDM
Vanilla
Time (ms/sample)
291.8
65.6
128.6
2770.5
361.3
65.2
Heat
Peak Memory (GB)
2,41
1.21
2.39
3.58
1.06
1.21
Time (ms/sample)
1371.0
780.6
211.2
8048.9
444.9
105.4
Burgers
Peak Memory (GB)
2.57
1.26
2.43
3.65
1.08
1.22
Time (ms/sample)
636.6
744.4
159.9
6578.7
474.7
80.7
Reaction-Diffusion
Peak Memory (GB)
2.97
1.72
2.93
4.33
1.43
1.43
Time (ms/sample)
350.7
3355.6
678.8
16447.2
338.7
343.8
Navier-Stokes
Peak Memory (GB)
7.39
2.03
3.90
5.85
1.93
2.02
Table 5: Runtime breakdown for PCFM showing the cost of sampling and the final projection step.
Dataset
Component
Time (ms/sample)
Percentage
Heat Equation
Sampling
286.6
98.22%
Final Projection
5.19
1.78%
Burgers Equation
Sampling
1358.8
99.1 99.11%%
Final Projection
12.2
0.89%
Reaction-Diffure Diffusion
Sampling
630.2
98.99%
Final Projection
6.40
1.01%
Navier-Stokes
Sampling
350.0
99.8%
Final Projection
0.7
0.2%
E.2 Tangent-Space Interpretation of the Projection Step
We provide additional theoretical context for the PCFM. In Proposition E. we formalize the projection step as tangent-space space update in Hilbert spaces and show that it corresponds an orthogonal projection onto the linearized constraint manifold, justifying its use for enforcing both linear and nonlinear constraints.
Proposition E. (Tangent-Space Projection Hilbert Spaces). Let H real Hilbert space and let h: 7 I be a Fréchet-differentaabl constraint operator Consider point un E H, and define the Jacobian = Dh(u1) L(H, )R,,, the bounded linear operator representing the Fréchet derivative.
20

---

Then, the update
Uproj - -JJ"
is the unique solution to the constrained minimization problem
min
u u1ll subject h(u1) J(u 0,
NEH
and corresponds to the orthogonal projection of un onto the affine subspace defined by the lineariza- tion of at un, i.e,, the tungent space to the constraint manifold u1
In the special case where h(u) Au - b is affine, with A E C(H, the update becomes the exact projection onto the feasible set {u H. Au b}.
Connection to ECI as a Special Case
The ECI (Extrapolate-Conm Coterlttnnterpllate)) framework [2] only considers and showcases empirical results constraints such fixed initial conditions, Dirichlet boundary values, and global mass integrals. These are all linear and non-overlapping, and are applied at isolated subsets of the solution field. To be used with the ECI framework, such constraints must admit closed-for oblique projections and act independently across disjoint regions. Let us denote any value or regional constraint [12] as a linear constraint huu) Au b, for some matrix E Rmxa vector b IR"". and for simplicity and comparison purposes, the ECI hyperparameters such as re-sampling is set to its default value (0) and mixing be 1.
In this regime, our Gauss-Newton update during PCFM guidance simplifies to an exact constraint satisfaction step for linear constraints:
Uproj = 41 - A A (AAT) - -),
which is precisely the "correction" " step ECI employs. Moreover, by disabling the relaxed penalty (A 0) and skipping intermediate corrections, PCFM reduces exactly ECI.
Hence, ECI special case of PCFM, tailored to simple linear constraints, while PCFM provides general and principled framework for enforcing arbitrary nonlinear equality constraints during generative sampling.
G Pre-trained Models with Functional Flow Matching
G.1 Flow Matching
39 formulates generative modeling as learning dimendeenn velocity field (u, T) that transports prior sample ug ~ to a target sample Ti. The induced flow or satisfies the ODE
Flow Matching (FM) 36
d
dr Φ5 (u) vo (1)) (u), (0(w)) tlo:
and defines a trajectory u(r) = Φ (uo) connecting To to T1. The learning objective minimizes the squared deviation between the model field and known target field é, evaluated along interpolants ur = T)uo + Tuj. This yields the standard Flow Matching loss:
(F) Br-2((01)) womeo, w1 [llva(u,, () (ur, 1]
(27)
where (u,, 7) is the true vector field, which is typically unknown.
Conditional Flow Matching (CFM). When the conditional velocity field ur ur(wi) := O tb, (uj) ) is
known in closed form, under optimal transport, one can minimize the conditional flow matching loss:
(CFM(() = Er (1,, 7) (w1 - uo))l
(28)
where u, = (1 - 7) + ru, is linear interpolant. This CFM objective forms the basis for most recent flow-based generative models due to its tractability and effectiveness.
21

---

G.2 Functional Flow Matching
To handle anfiite-diimeniional generative tasks such as PDE solutions. Functional Flow Matching (FFM) [[] extends CFM to Hilbert spaces U of functions u : X' + R. The FFM flow th, satisfies
d
dr Vr (w) - uo (4), ug μo:
where μo is a base measure function space. Under regularity assumptions, FFM minimizes an analogous loss over the path of measures:
(0) Er.no.w1 [Ivo(urs 7) (wi - mo) =0))
(29)
where u, (1 r) + Tuı and the L -orm is used to evaluate function-valued velocities over the domain X. This loss generalizes CFM to the functional setting, and serves as the pretraining objective for PCFM in this work.
H Sampling Setups for PCFM and Other Baselines
For all comparisons, we adopt the explicit Euler integration scheme unless otherwise specified. We use 100 Euler update (flow matching) steps for heat and Navier-Stokes, and 200 steps for Reaction- Diffusion and Burgers (IC or BC).
Vanilla FFM. We perform unconstrained generation by integrating the learned flow vector field using explicit Euler steps, following the standard flow matching procedure 4 IT. At each step, the model applies the forward update U++
dt + At (t, up) without any constraint enforcement.
PCFM (Ours). PCFM performs explicit Euler integration combined with constraint correction every step. We apply Newton update to project intermediate states onto the constraint manifold. Additionally, optionally apply guided interpolation with À 1.0, step size 0.01, and 20 refinement steps to further refine interpolated states, although we use À = O by default (no interpolation guidance) for computational efficiency (see Appendix L.3 for further ablations
A). Interpolation guidance is only applied for the heat equation case to obtain results in Table For Navier-Stokes, we follow the stochastic interpolant idea and randomize noise batches 37 (i.e., different noise samples over batches). Finally, after all flow matching steps, we apply a final projection on the solution to enforce the constraints.
ECL. We follow the ECI sampling procedure introduced by Cheng et al. [12] where it performs iterative
steps well as several mixing and noise resampling steps throughout the trajectory. For the simpler heat equation case, we did not do mixing and noise resampling. For more challenging PDEs Navier-Stokes, Reaction-Diffume
and Burgers), we apply Amix = 5 mixing updates per Euler step, and resample the Gaussian prior every 5 steps. We adopt their value enforcements in IC and Dirichlet BC cases, and also constant mass integral enforcement in the heat and Navier-Stokes problems, where the periodic boundary conditions lead to trivial linear constraints.
D-Flow. Following Ben-Hamu et al. [35], D-Flow differentiates through an unrolled Euler integra- tion process. Following their paper and the set-up in Cheng et al. [1), we optimize the initial noise via LBFGS with 20 iterations and learning rate of 1.0. The optimization minimizes the constraint violation at the final state, requiring gradient computation through the full unrolled trajectory. We adopt the adjoint method shipped from torchdiffeg (38) to differentiate through the ODE solver. Depending the problem,
adopt an IC BC loss well as PINN-loss based the differential form of the PDE to form the constraint loss function to be optimized. We use 10 iterations and a learning rate of 10 for the Navier-Stoker dataset to avoid NaNs in the optimization loop.
DiffusionPDE. We adopt a gradient-guil -uided sampling method proposed by Huang al. 10 on the same pretrained FFM model. The generation process is augmented with explicit constraint correction where a composite loss of IC/BC and PINN [22] (differential form of the PDE), like the D-Flow set-up, is used to guide the vector field at each flow step. We apply correction coefficient r/ 1.0 and set the PINN loss weight to 10-2 2 for stable generation. We find that a higher PINN weight or higher learning rate lead to unstable generation.
22

---

PDM. We implement the PDM algorithm proposed by Christopher et al. [49 for constrained generation. Although PDM also employs projection during sampling, its assumptions and modeling framework differ fundamentally from those of PCFM. For a fair comparison, we additionally im- plemented PDM-style ablation within our flow models using the proposed projection step. While PDM achieves comparable performance tently outperforms on more complex nonlinear dynamics, suggesting that PDM's strict per-step projections over-constrain the sampling trajectory and degrade generation quality.
simpler systems such as the Heat equation, PCFM consis-
I Batched Differentiable Solver for Nonlinear Constraints
We describe solver for projecting predicted sample un IR" onto the nonlinear constraint manifold M == {u E IR": h(u) = 0}, where R + IR".. The projection formulated as the
constrained optimization problem:
min H| - un 11 subject to h(u) = O.
(30)
wERn
The corresponding Lagrangian is: [(u, = g1lue 21 + AATh(u),
with first-order optimality conditions:
VuL(u, A) = u un + J(u) > ==
h(u)= = 0,
where J(u): == Vh(u) E Rmx" is the constraint Jacobian.
These yield the full nonlinear KKT system:
ui + + 0, h(uu= =.
(31)
Newton-Based Update. At iteration k we linearize the KKT system around u (*),, ( and solve for updates 8u, 5X. The full Newton step solves:
{ (h((((() H(lk))
)
6u
(k)
+ +(((k))
(32)
82
h(u(k))
Here, the upper-left block contains the Hessian of the Lagrangian:
m
V2 C() (()) = I
(0) (((((())))
After solving, we update the primal and dual variables:
n (k+1) = u (8) + Su, (M+11 (4) + 5X.
This Newton system converges quadratically under standard regularity assumptions (e.g.. full-rank Jacobian and Lipschitz-comtinnous second derivatives). In practice, we often omit the second-order terms to yield Gauss-Newton approximation that is more stable and efficient in high-dimensional settings.
Approximate KKT Solve via Schur Complement. For inference-ting projection, we adopt simplified and batched update using the Schur complement [69]. At each iteration, we set 0, and solve only for the primal update. Eliminating du from Equation (32 we obtain:
(JJT) X A= -h(u,
du = J J
(33)
This gives the Neesso-ewton-stee style update:
u + - J (JJJ)-1 (h(u) + [[1 - u))
(34)
We iterate this procedure until convergence or until the constraint residual [[h(w]]] falls below set tolerance. The matrix JJJ E mm is small and typically well-conditioned for local or sparse constraints, enabling efficient solves.
Restarting the dual variable with X 0 at each iteration avoids stale gradient accumulation and improves numerical stability. This leads to a robust and memory-efficient projection routine that supports batched execution and reverse-mode differenniaaion..
23

---

Batched and Differentiable Implementation. We implement the solver in batched and differ- entiable fashion support inference across samples. For a batch of inputs {ui] we evaluate Jacobians J', residuals h(w"), and solve the corresponding Schur systems in parallel using vector- ized operations and autodiff-compatible backends (e.g., PyTorch with batched Cholesky or linear solvers) 71. Potentially, specific GPU kernels can be built for ODE integration to ensure optimal performance 38 72].
Computational Complexity. The per-sample cost includes:
(m)") for computing J and J
(m)) for solving the Schur complement system,
(n) for applying the update.
Since typically m < n. the overall cost scales as (n)) per sample.
J Evaluation Metrics
We evaluate each method using 512 generated samples for all ID problems and 100 samples for the 2D Navier-Stokes problem. For each setting, solutions with a fixed or BC configuration used during sampling for direct comparison.
an equivalent number of ground truth PDE
MMSE and SMSE. Following Kerrigan et al. T. Cheng et al. [2]. we evaluate generation fidelity with mean of the MSE and the standard deviation of the MSE as:
MMSE
IlMgen Mgg
SMSE
||0
0gil
where Mgens Mg: denote the mean across generated and true PDE solutions, and gens Qgt the standard deviations.
To evaluate physical consistency, compute the le norm of residuals from constraint functions R, applied to each sample @(m), then average across all N samples:
Constraint Error.
CE(+)
Σ R.
BC, CL).
where the residuals are computed following Appendix violation.
D to measure the constraint or conservation
Fréchet Poseidon Distance (FPD). To assess distributional similarity beyond mean and variance, we adopt the Fréchet Poseidon Distance (FPD) introduced in Cheng et al. (12), where we measure the Fréchet distance between the hidden state distributions extracted from pretrained foundation model (Poseidon (67] applied to both generated and true solutions. We pass both generated and true solutions through the Poseidon base model (157M parameters) and extract the last hidden activations of the encoder to eventually obtain 784-dimension feature vector for FPD calculation:
FPD" IIl =all + Tr (3) 2
(35)
where μ1,, Σ₁ and μ12 Σg are the empirical mean and covariance of the Poseidon embeddings from the generated and true solutions' distributions, respectively. FPD is computed either per frame u(x v) for 2D PDEs and averaged over all frames (across () or over the full spatiotemporal solution u(æ, L) for ID PDEs.
24

---

K Further Results
Deattva
Generated solution mean
Generated solution std dev
residual mean
(00-1)
PCNN
DlustonPO
PON
residual mean
Figure 5: Solution profiles for the Inviscid Burgers equation with fixed BC. We plot the various constraint guidance methods and compare the mean solution profile and standard deviation. While PCFM yields slightly worse MMSE and SMSE and better FPD, ensures global mass conservation and maintains low constraint errors for both Dirichlet and Neumann BCs over time.
DifusioPPOE
OF
FFM
PDM
Generated solution mean
Generated solution std dev
64
Figure 6: Solution profiles for the Heat equation with fixed IC. We plot the various constraint guidance methods and compare the mean solution profile and standard deviation. PCFM outperforms all other methods by visually being the most similar to the ground truth.
25

---

MH
Drhabon05
Generated
solution
Ground Truth
PCFM
FFM
Dflow
EC
DiflusionPDE
PDM
250
2.35
solution
std dev
Figure 7: Alternative view of the Reaction-Diifussiin equation with fixed IC. We plot the pointwise mean and standard deviation across generated samples for each method. PCFM outperforms all other methods by visually being the most similar to the ground truth.
Generated solution mean
Generated solution std dev
Figure 8: Alternative view of the Inviscid Burgers equation with fixed IC. We plot the pointwise mean and standard deviation across generated samples for each method. PCFM outperforms all other methods by visually being the most similar to the ground truth.
L Further Ablations
L.1 Effect of Intermediate Projections
Physical constraints such as mass conservation or initial conditions are typically imposed only on the final state. However, find that applying projections at intermediate steps improves flow consistency and numerical stability. This parallels practices in diffusion and flow-matching models, where intermediate states, often initialized noise, need not satisfy constraints carry semantic meaning.
To assess this, we perform ablation presented in Table in which projections applied only at final step, omitting intermediate corrections. observe that trajectories drift farther from the constraint manifold, resulting larger corrective updates, slower convergence, and lower accuracy.
Intermediate projections, by contrast, act as continuous guidance, progressively steering samples toward the feasible manifold while preserving diversity. This avoids over-constrainin noisy interme- diates and yields more accurate and stable final solutions.
L.2 Total Variation (TV) Constraints the Heat Equation
Total Variation Diminishing (TVD) constraints encode the principle that diffusive systems, such as those governed by the heat equation, smooth spatial fluctuations time. Mathematically. this is expressed by the fact that the total variation of the solution should not increase in time:
du dr
TV(u(())
dx is non-increasing in f.
To enforce this property in a data-driven or generative setting, we adopt a hard constraint that relates the total variation at final time T' to the total variation the initial time 0. Specifically,
26

---

Table Ablation study on intermediate projections. Applying projections only at the final step leads to larger constraint deviations and higher errors, while PCFM's intermediate corrections improve stability and accuracy across all datasets.
Dataset
Method
MMSE/ 10-2
SMSE/ 10
CE (IC)
CE (CL)
Heat
Project final only PCFM
4.31
3.27
0.241
0.937
0
0
Burgers IC
Project final only PCFM
10.59
6.90
0
4.05
0.052
0.272
Reaction-Diffusion
Project final only
2.82
2,48
PCFM
0.026
0.583
Navier-Stokes
Project final only
16.47
8.60
PCFM
4.59
4.17
0
impose the condition:
TV(ur) 7 -V-TVuo)),
where (0, 1) is decay factor that can either be fixed (e.g.. 0.5) or estimated from unprojected samples. This constraint can implemented using discrete spatial differences as:
TV(u,)
~
1441,
applied at selected time slices (typically = 0 and = 7). We encode the constraint differentiable function:
(T) == TV(wr) TV(wo) = 0.
This constraint is compatible with conservation and initial condition enforcement, and can be combined with other physically informed conditions such as energy decay or curvature regularity. In practice, observe that enforcing TVD constraints reduces spurious oscillations in the generative output and improves structural fidelity without significantly altering statistical accuracy.
Power MSE.
To further quantify the effect of this constraint, we introduce the Power MSE metric, which measures the deviation in spatial frequency content between the generated and ground truth solutions. Formally, for a model's 's mean output ū(x t), we compute the spatial power spectrum [ū(k; 1)(² and average over time to obtain a per-frequency energy profile. The Power MSE is then defined as:
The Σ 2(1gen (k) -
Power MSE
This metric captures discrepancies in frequency modes and is especially sensitive to unphysical sharpness or over-smoothing. Variants of power-spectrum-basee error metries have been used in turbulence modeling (73), climate downscaling [74), and Fourier neural operators training [23]. In our setting, we find that Power MSE effectively complements pointwise metries (e.g.. MMSE, CE), as it evaluates structural fidelity in the spectral domain. As shown in Table 7 this metric reveals how enforcing TVD constraint not only improves calibration but also leads closer alignment in frequency space.
Effect of TVD Constraint. To isolate the benefit of total variation control, mirror the ECI setting by setting the constraint penalty weight to zero and enforcing only initial condition and mass conservation constraints, When we additionally include the TVD constraint (y 0.3) in our PCFM model (with no penalty- based regularization), we observe consistent improvement in structural fidelity. This is most evident in the Power MSE metric, where PCFM achieves the lowest error among all baselines (488.42 vs. 500.12 for ECI), indicating better alignment with the true spatial energy distribution. We also see gains MMSE and SMSE, without introducing constraint violations (e.g.. CE remains zero). This confirms that the TVDD constraint acts as meaningful structural prior even in the absence of learned penalties, improving generative realism without sacrificing calibration or efficiency.
27

---

Table 7: Test metrics on the Heat Equation dataset with TVD constraints. Lower is better for all metrics.
Metric
PCFM
ECI
DiffusionPDE
D-Flow
FFM
MMSE/ / 10
0.684
0.697
4.49
1.97
4.56
SMSE / /10-2
0.962
0.973
3,93
1.14
3.51
CE (IC) / 10-2
0
0
599
102
579
CE (CL) 10-2
0
0
2.06
64.8
2.11
FPD
1.28
1.34
1.70
2.70
1.77
Power MSE
488.42
500.12
1420.92
3996.11
932.20
L.3 Effect of Relaxed Constraint Correction with Flow Matching Steps
MMSE SMSE
SMSE
0.010
Eroor
Emoor
MMSE
SMAE
0.004
6.000
0.002
Ponalty weight
Peaalty weight 7
Figure 9: Effect of penalty weight À MMSE and SMSE for the Reaction-Diffusion dataset. Left: 10 flow matching steps. Right: 100 flow matching steps.
Figure 9 shows the effect of the relaxed constraint penalty weight À on the generation error, mea- sured via MMSE and SMSE, for the Reaction-Diffusee dataset. We evaluate on 128 100 spatial-temporal grid, described in Appendix B.3 and solve the relaxed optimization objective in Equation (6) using the Adam optimizer.
With only 10 flow steps (left), increasing À significantly improves performance, as the relaxed constraint correction effectively compensates for the coarser integration of the flow. In contrast, with 100 flow steps (right). the flow is more precise, and additional penalty yields only marginal benefit, This illustrates that relaxed correction particularly valuable when inference is constrained a small number of steps for instance, in scenarios where evaluating the vector field vo is computationally expensive.
However, overly large values of À can harm performance by distorting the reverse update and breaking alignment with the OT interpolant (see Proposition as discussed in Equation (5). Hence, choosing requires balancing constraint enforcement with consistency along the learned generative path.
L.4 Evaluation on Generating Dutts--Distributtin Samples
A key advantage of inference-time constraint enforcement is its ability to generalize beyond the constraints implicitly present the training data. To assess this capability, we evaluate PCFM under dustof-distribuutio (OOD) constraint settings, where the pretrained model encounters constraint realizations that were never observed during training. Specifically, we design an OOD experiment based the Burgers equation, focusing shifts in the initial condition (IC). The pretrained flow- matching model is trained on standard dataset in which the initial sigmoid profile
u(x, 0) exp(-a(x= - mo)]
centered at locations 20 E [0.4, 0.6]. For testing, we construct a new OOD dataset where oy sampled from the non-overlapping intervals zo e [0.0, 0.2] and E [0.8, 1.0). This shift introduces
28

---

Table 8: Performance out-of-distribution (OOD) dataset for the Burgers equation. PCFM achieves the lowest errors while exactly satisfying both mnttial-coniitinn (IC) and conservation-lav (CL) constraints, demonstrating robust generalization beyond maaning-distribune constraints.
Method
MMSE/ 10
SMSE 10-2
CE (IC)
CE (CL)
PCFM
1.83
1.40
0
ECI
13.30
5.65
3.64
D-Flow
15.59
7.82
6.00
0.092
Vanilla
18.84
7.88
6.37
0.069
DiffusionPDE
19.49
8.03
6.42
0.062
PDM
206.6
2.81
0
initial conditions that lie outside the training distribution, thereby evaluating each method's ability to adapt to unseen constraint configurations.
All models are evaluated retraining. We employ 200 flow-matching steps for the Burgers problem to ensure stable integration and consistent comparison across methods. Each approach enforces the same IC and conservation constraints, providing a fair assessment of generalization ability.
zero-shot setting, reusing the pretrained FFM backbone without
The results, summarized in Table show that PCFM achieves the lowest mean and sample-wise mean squared errors (MMSE and SMSE) while exactly satisfying both IC and conservation-law constraints, under this challenging OOD regime. In contrast, unconstrained partially constrained baselines exhibit significant violations, confirming that inference-time projection enables robust generalization unseen constraint configurations.
29

---

NeurIPS Paper Checklist
The checklist is designed to encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No], or [NA].
[NA] means either that the question is Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justifiicton right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors in their evaluation. While "[Yes] is generally preferable "[No] ", it is perfectly acceptable answer "[No] provided a proper justification is given (e.g., "error bars are not reported because it would be too computationally expensive" "we were unable find the license for the dataset we used"). In general, answering "[No] or "[NA] is not grounds for rejection. While the questions are phrased binary way. we acknowledge that the answer is often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendíx. If you answer [Yes] a question, in the justification please point to the section(s) where related material for the question can be found.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The authors study the problem of imposing arbitrary nonlinear hard constraints in flow-based generative models a zero-shot manner. The authors list their contributions appropriately in the abstract and introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should elearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
The authors list the limitations of their work in the Conclusion section.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
30

---

The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, c.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-texl system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full of assumptions and a complete (and correct) proof?
Answer: [NA]
Justification:
The paper does not include any major theoretical result.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: The authors provide the details of their baselines (Appendix H) test setup
(Appendix B Appendix C and evaluation metrics (Appendix[].
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data provided or not.
31

---

If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g.. to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The paper provides access to the data and code, with details of reproducibility in Appendix
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cm cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
//nips ltaips.cc/public/guuido htt.ggudddde/oeeuubiii.i..ee
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
32

---

Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: The whole training and test details are described in Appendix H Appendix []] and Appendix !
Appendix
B
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We report metrics such as Standard Mean Square Error (SMSE) and variance in the mass residual mean plots in Figure 2 and Figure
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The compute resources are listed in Appendix C
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
33

---

The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: The authors confirm that the research conducted in the paper conforms, in every respect, with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e-g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We discuss the broader impact of our work in the conclusion section.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deploymeet of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied partieular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out. For example, it is legitimate point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML..
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses such risks.
Guidelines:
The answer NA means that the paper poses such risks.
34

---

Released models that have high risk for misuse or dual-use use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The authors cite the original paper that produced the code package dataset with proper attribution.
Guidelines:
The answer NA means that the paper does existing
The authors should cite the original paper that produced the code package dataset.
The authors should state which version asset used and, possible, include URL.
The name of the license (c.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g.. website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, batasets.paaerssttt paperswi thcode. cow/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release any assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects.
35

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The research mentioned in the paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMS if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: We only use LLMs for writing, editing, or formatting purposes.
Guidelines:
The answer NA means that the core method development in this research does involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: /tte::///rurrre iis cc/ /cc/Confeeeencee2225////////// for what should or should be described.
36


Reviews:
[{"review_id": "DUXkxJYHpT", "reviewer": "Reviewer_8V7Y", "summary": "This work proposes Physics-Constrained Flow Matching (PCFM), a zero-shot inference framework for enforcing arbitrary nonlinear equality constraints in pretrained flow-based generative models. PCFM combines flow-based integration, tangent-space projection, and relaxed penalty corrections to strictly enforce constraints while staying consistent with the learned generative trajectory. This approach supports both local and global constraints, including conservation laws and consistency constraints, without requiring retraining or architectural changes. Empirically, PCFM outperforms state-of-the-art baselines across diverse PDEs, including systems with shocks and nonlinear dynamics, achieving lower error and exact constraint satisfaction.", "questions": "1.  A practical comparison of wall-clock inference times against the baselines would provide a more concrete sense of the computational efficiency of PCFM.", "limitations": "N/A", "rating": 5, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "I do not have enough background to review the paper. Here are a few comments:\n\n**Strengths**\n\n1. The \"shoot-project-reverse\" cycle designed in this work looks novel. This work presents an effective solution for satisfying the strict physical laws required when applying generative models to physical systems. It enforces complex, nonlinear equality constraints as hard requirements in a general way during the inference stage alone, without requiring to retrain the model or modify its architecture to introduce biases.\n\n2. Empirical results in Table 3 looks promising. The proposed method seems much better than the existing counterparts. \n\n**Weaknesses**\n\n1.  In the step of forward shooting and projection, the cost depends on forming the constraint Jacobian, $J = \\nabla h(u_1)$. For high-dimensional, complex, and non-local physical constraints, does constructing the Jacobian $J$ lead to a computational explosion?", "was_revised": true, "final_justification": "Thanks for the authors' rebuttal. My questions are clearly answered. I think the paper is technically solid with significant contribution. I keep the final rating unchanged.", "rebuttal": "We thank the reviewer for their thoughtful assessment and for highlighting both the novelty of our shoot–project–reverse cycle and the empirical strength of our method across challenging physical domains. We are encouraged that the proposed framework was found effective and promising for enforcing strict physical constraints at inference time.\n\n## R1.1: Computational Cost of Jacobian Construction\n\nWe appreciate the reviewer raising this important question regarding scalability.\nIn practice, the constraint Jacobian Jh(x) is not constructed explicitly. Instead, we compute it using batched reverse-mode automatic differentiation, enabled via vectorized mapping (e.g., vmap in PyTorch). This allows us to efficiently compute the full Jacobian across multiple constraint functions in parallel. Since each row of the Jacobian corresponds to a reverse-mode gradient computation, this formulation is highly parallelizable on GPUs and scales well even when the state dimension is large.\n\nFor solving the projection step, employing iterative solvers (e.g., conjugate gradient, GMRES) on the Gauss–Newton system will require only Jacobian-vector products (Jv) rather than the full matrix. These products are efficiently handled by automatic differentiation and avoid the memory costs associated with storing large dense Jacobians.\n\nAdditionally, our approach is compatible with Jacobian-Free Newton–Krylov (JFNK) methods, where even the Jv products are approximated via finite-difference evaluations. This further reduces memory footprint and makes the solver suitable for extremely high-dimensional or non-local constraint settings.\n\nIn summary, the projection step does not pose a computational bottleneck, and we find it to be practical and scalable in our experiments. \n\n## R1.2: Wall-clock inference times\n\nWe provide additional results, measuring our inference times compared with the existing baselines. To compare fairly on time and peak memory, we use the same batch size for all methods and loop through the test set to record the per sample time and peak memory. We use 100 flow matching steps across all methods and datasets for benchmarking here. We use a single Nvidia V100 GPU with 32GB ram for all experiments. Batch sizes are 32 for heat, reaction diffusion, and Burgers problems and 8 for Navier-Stokes. We note that DFlow runs out of memory even for a batch size of 1 on the V100 chip, thus we improved the algorithm by including an adjoint ODE solve. For PDM, we adopt our constraint projection to every intermediate state along the flow matching trajectory. We report the average per-sample inference time and peak memory usage. \n\nPCFM incurs a modest runtime overhead compared to unconstrained baselines like vanilla FFM, due to projections to satisfy nonlinear hard constraints exactly. However, it achieves 4-6x speedups over D-Flow, which optimizes noise by backpropagating through the full ODE unroll - which is highly memory-intensive and slow. For Heat and Navier-Stokes datasets, the constraints are linear with well-defined structure, the projections converge in a single step, evident also in PDM's runtime. DiffusionPDE is faster than PCFM but it compromises on generative quality and struggles to enforce hard constraint, particularly for nonlinear problems like Burgers and Reaction-Diffusion. Similarly, ECI's mixing iterations can increase runtime overhead when tuned for more complex problems but it fails to achieve hard constraint satisfation for nonlinear problems.\n\nProfiling in Table B reveals 96-99% of PCFM's time is in intermediate projection steps, with final projection at 1-3% of the total runtime. \n\n\n### Table A: Runtime and Memory Comparison\n\n| **Dataset**          | **Metric**       | **PCFM** | **ECI** | **DiffusionPDE** | **D-Flow** |  **PDM**  | **Vanilla** |\n| -------------------- | ---------------- | :------: | :-----: | :--------------: | :--------: | :-------: | :---------: |\n| **Heat**    | Time (ms/sample) |  291.8  |  65.6  |      128.6      |   2770.5  |   361.3  |  65.2  |\n|                      | Peak Memory (GB) |   2.41   |   1.21  |       2.39       |    3.58    |    1.06   |   1.21  |\n| **Burgers** | Time (ms/sample) |  1371.0 |  780.6 |      211.2      |   8048.9  |   444.9  |  105.4 |\n|                      | Peak Memory (GB) |   2.57   |   1.26  |       2.43       |    3.65    |    1.08   |   1.22  |\n| **Reaction-Diffusion** | Time (ms/sample) | 636.6 | 744.4 | 159.9     |  6578.7 |   474.7  |  80.7 |\n|                      | Peak Memory (GB) |   2.97  |  1.72  |     2.93       |   4.33    |    1.43 |   1.43  |\n| **Navier–Stokes**    | Time (ms/sample) | 350.7 | 3355.6 | 678.8      | 16447.2   | 338.7  | 343.8 |\n|     | Peak Memory (GB) | 7.39  |   2.03    |  3.90     |  5.85  |  1.93   | 2.02 |  \n---\n\n\n### Table B: PCFM Runtime Breakdown\n\n\n| **Dataset**          | **Component**                           | **Time (ms/sample)** | **Percentage** |\n| -------------------- | --------------------------------------- | ------------- | -------------- |\n| **Heat Equation**    | Sampling  | 286.6        | 98.22%         |\n|                      | Final Projection        | 5.19          | 1.78%          |\n|                      | **Total**                               | **291.8**    | 100.00%   |\n| **Burgers Equation** | Sampling | 1358.8        | 99.11%         |\n|                      | Final Projection         | 12.2          | 0.89%          |\n|                      | **Total**                               | **1371.0**    | 100.00%    |\n| **Reaction-Diffusion**    | Sampling |   630.2   |   98.99%       |\n|     | Final Projection   |    6.40     |       1.01%   |\n|        | **Total**                       | **636.6** | 100.00%    |\n| **Navier–Stokes**    | Sampling |    350.0    |     99.8%    |\n|                      | Final Projection    |  0.7   |     0.2%       |\n|                      | **Total**            |  **350.7**   |  100.00%    |"}, {"review_id": "MD52hJmku8", "reviewer": "Reviewer_aA1E", "summary": "This paper proposes Physics-Constrained Flow Matching (PCFM), a zero-shot inference framework to enforce physical constraints  in flow-based generative models. PCFM corrects intermediate and final samples by projecting onto the constraint manifold, without requiring retraining or access to constraint gradients during training.", "questions": "See Weaknesses", "limitations": "See Weaknesses", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n- The shock fronts phenomenon when solving the Burgers's equation in Fig.1 are quite impressive.\n- The paper identifies a core challenge in generative surrogate modeling for physics\n- The writing is clear and easy to follow.\n\nWeaknesses:\n\n- I believe the term *hard constraint* may be an overstatement. In my view, a constraint can only be classified as *hard* if it is mathematically satisfied exactly (i.e., equals zero). Since your method merely reduces the violation to a small value rather than eliminating it entirely, it might be more appropriate to refer to it as a *soft constraint* or an *approximate hard constraint*.\n- I am curious about the necessity of projecting the guess from intermediate flow states onto the constraint manifold. Given that physical constraints are typically imposed only on the final state, it is unclear what benefits such intermediate projections provide. Moreover, unless the constraint function h is affine, your projection step does not generally yield an exact solution on the target physics manifold. A more detailed justification would be helpful.\n- In line 54, you mention that some methods require backpropagation through expensive PDE operators. However, your projection step also appears to be computationally intensive. It would be beneficial to include a comparison of computational time and memory consumption between your method and those baselines, to better contextualize the trade-offs involved.\n- Could you provide additional qualitative comparisons between your method and all baseline models on the Navier-Stokes dataset? The other datasets used in your paper are relatively simple, which may limit the generalizability of your conclusions. Furthermore, experiments on higher-resolution or more complex datasets are highly recommended to more robustly demonstrate the effectiveness of your approach.", "was_revised": true, "final_justification": "I recommend this paper to accept when the authors promise to avoid using term \"hard constraints\" incorrectly.", "rebuttal": "We thank the reviewer for their constructive feedback and for appreciating the clarity of the writing, the motivation behind enforcing physical constraints in generative models, and the strength of the shock front results in Figure 1. Below, we address each of the key concerns in detail.\n\n## R2-1: Use of the Term “Hard Constraint”\nWe appreciate the reviewer’s attention to terminological precision. In our work, we enforce equality constraints via Gauss–Newton projection, with residual norms consistently reduced to within machine epsilon. In Float32, this corresponds to violations on the order of $10^{-7}$, and with Float64, the residuals drop to $\\sim 10^{-15}$. This level of constraint satisfaction aligns with the definition of hard constraints in the numerical optimization literature, where constraints are considered “hard” if they are satisfied up to solver and machine precision [1].\nWe will clarify this usage in the revised manuscript and explicitly state that our method achieves hard constraint satisfaction to the limits of numerical precision.\n    \n## R2-2: Necessity and Role of Intermediate Projections\n\nWhile constraints in physical systems (e.g., mass conservation or initial conditions) are often imposed only on the final state, we find that intermediate projections offer clear benefits for flow consistency and numerical stability. This aligns with practices in diffusion and flow matching literature, where intermediate states—often initialized as noise—are not expected to satisfy constraints or carry semantic meaning (e.g., intermediate samples in Stable Diffusion are noisy and not valid images).\n\nTo assess this, we conducted an ablation in which projections were applied only at the final step, omitting intermediate corrections. Using the same setup (100 flow steps for Heat; 200 for Burgers, RD, and NS), we observed that trajectories drift farther from the constraint manifold, leading to larger corrective steps, slower convergence, and reduced accuracy. This confirms that relying solely on final correction is suboptimal, especially for nonlinear or global constraints.\n\nIn contrast, intermediate projections act as **lightweight continuous guidance**, gradually steering samples toward the feasible set while preserving diversity. This avoids overly rigid enforcement on noisy intermediate states and improves final performance. Our results validate this approach, and we will include the ablation and clarify the motivation in the revision.\n\n### Table A: Comparisons between only projecting final step and PCFM\n| Dataset | Method    | MMSE      | SMSE      | CE (IC) | CE (CL) |\n|---------|-----------|-----------|-----------|--------------|---------------|\n| Heat    | project final only  | 0.0431 |  0.0327 | 0 |  0 |\n| Heat    | PCFM      |    **0.00241**    |   **0.00937**       |          0    |        0       |\n| Burgers | project final only  | 0.1059 | 0.0690 | 0 | 4.0515 |\n| Burgers | PCFM      |   **0.00052**        |     **0.00272**      |      0      |      0         |\n| Reaction-Diffusion      | project final only  | 0.0282 | 0.0248 | 0 | 0 |\n| Reaction-Diffusion      | PCFM      |    **0.00026**       |     **0.00583**      |     0       |        0       |\n| Navier-Stokes      | project final only  | 0.1647 | 0.0860  | 0 | 0  |\n| Navier-Stokes      | PCFM      |    **0.0459**       |   **0.0417**       |    0          |      0         |\n\n## R2-3: Computational Cost of Projection\n\nThe reviewer correctly notes that constraint projections can introduce computational overhead. In PCFM, the projection step solves a Gauss–Newton system using iterative solvers (e.g., CG, GMRES) that require only Jacobian-vector products (Jv), computed efficiently via reverse-mode autodiff. We further accelerate this using batched vectorized maps (e.g., `vmap` in PyTorch) for GPU parallelism.\n\nOur implementation also supports **Jacobian-Free Newton–Krylov (JFNK)** methods, which approximate directional derivatives via finite differences, eliminating the need for explicit Jacobian construction and scaling well to high-dimensional and nonlocal constraints.\n\nFor fair benchmarking, we use consistent batch sizes and 100 flow steps across all methods and datasets, evaluating average per-sample runtime and peak memory on a single NVIDIA V100 GPU. DFlow, which backpropagates through the entire ODE solve, is significantly slower and more memory-intensive; even with adjoint optimization, PCFM is 4–6× faster on average. In contrast, unconstrained baselines like vanilla FFM are faster but cannot guarantee constraint satisfaction. DiffusionPDE is faster than PCFM but sacrifices constraint enforcement and quality on nonlinear tasks. ECI’s runtime varies with mixing steps but fails to enforce nonlinear constraints robustly.\n\nProfiling (Table C) shows the final projection step contributes only 1–3% of total time; most cost is due to sampling itself. Overall, PCFM introduces a modest and practical overhead that is necessary to enforce hard nonlinear constraints, and we will summarize this analysis in the main paper.\n\n### Table B: Runtime and Memory Comparison\n\n| **Dataset**          | **Metric**       | **PCFM** | **ECI** | **DiffusionPDE** | **D-Flow** |  **PDM**  | **Vanilla** |\n| -------------------- | ---------------- | :------: | :-----: | :--------------: | :--------: | :-------: | :---------: |\n| **Heat**    | Time (ms/sample) |  291.8  |  65.6  |      128.6      |   2770.5  |   361.3  |  65.2  |\n|                      | Peak Memory (GB) |   2.41   |   1.21  |       2.39       |    3.58    |    1.06   |   1.21  |\n| **Burgers** | Time (ms/sample) |  1371.0 |  780.6 |      211.2      |   8048.9  |   444.9  |  105.4 |\n|                      | Peak Memory (GB) |   2.57   |   1.26  |       2.43       |    3.65    |    1.08   |   1.22  |\n| **Reaction-Diffusion** | Time (ms/sample) | 636.6 | 744.4 | 159.9     |  6578.7 |   474.7  |  80.7 |\n|                      | Peak Memory (GB) |   2.97  |  1.72  |     2.93       |   4.33    |    1.43 |   1.43  |\n| **Navier–Stokes**    | Time (ms/sample) | 350.7 | 3355.6 | 678.8      | 16447.2   | 338.7  | 343.8 |\n|     | Peak Memory (GB) | 7.39  |   2.03    |  3.90     |  5.85  |  1.93   | 2.02 |  \n\n### Table C: PCFM Runtime Breakdown\n\n| **Dataset**          | **Component**                           | **Time (ms/sample)** | **Percentage** |\n| -------------------- | --------------------------------------- | ------------- | -------------- |\n| **Heat Equation**    | Sampling  | 286.6        | 98.22%         |\n|                      | Final Projection        | 5.19          | 1.78%          |\n|                      | **Total**                               | **291.8**    | 100.00%   |\n| **Burgers Equation** | Sampling | 1358.8        | 99.11%         |\n|                      | Final Projection         | 12.2          | 0.89%          |\n|                      | **Total**                               | **1371.0**    | 100.00%    |\n| **Reaction-Diffusion**    | Sampling |   630.2   |   98.99%       |\n|     | Final Projection   |    6.40     |       1.01%   |\n|        | **Total**                       | **636.6** | 100.00%    |\n| **Navier–Stokes**    | Sampling |    350.0    |     99.8%    |\n|                      | Final Projection    |  0.7   |     0.2%       |\n|                      | **Total**            |  **350.7**   |  100.00%    |\n\n\n## R2-4: Generalizability and Experimental Scope\n    \nWe agree that demonstrating applicability beyond PDE-based systems would broaden the paper’s scope. Our primary focus is on scientific generative modeling, where physical constraints play a central role and quantitative evaluations are well-defined. To ensure comparability with prior work, we focused on **standard benchmark PDE datasets**, including **smooth equations** such as Heat and Reaction–Diffusion, and **nonlinear, shock-forming systems** such as Burgers and Navier–Stokes. These settings are commonly used in the constrained generative modeling literature \\[2, 3, 4], neural operator learning \\[5, 6], and foundation PDE models \\[7, 8]. In addition, we also provide runtimes to all our datasets and add an additional baseline (PDM) \\[9\\], including on the three-dimensional Navier-Stokes dataset as well. The underlying FNO architecture in the PCFM can also be naturally adopted zero-shot super-resolution \\[5\\]. Nevertheless, we recognize the value of extending to more complex or real-world settings.\n    \n[1] Nocedal, Jorge, and Stephen J. Wright. Numerical optimization. New York, NY: Springer New York, 2006.\n    \n[2] Cheng, Chaoran, et al. \"Gradient-free generation for hard-constrained systems.\" (2025).\n    \n[3] Huang, Jiahe, et al. *DiffusionPDE: Generative PDE-solving under partial observation.* arXiv:2406.17763, 2024.\n\n[4] Ben-Hamu, Heli, et al. *D-flow: Differentiating through flows for controlled generation.* arXiv:2402.14017, 2024.\n\n[5] Li, Zongyi, et al. *Fourier neural operator for parametric partial differential equations.* arXiv:2010.08895, 2020.\n\n[6] Li, Zongyi, et al. *Physics-informed neural operator for learning partial differential equations.* arXiv:2111.03794, 2021.\n\n[7] Herde, Maximilian, et al. *Poseidon: Efficient Foundation Models for PDEs.* arXiv:2405.19101, 2024.\n\n[8] Sun, Jingmin, et al. *Towards a Foundation Model for Partial Differential Equation: Multi-Operator Learning and Extrapolation.* arXiv:2404.12355, 2024.\n\n[9] Christopher, Jacob K., Stephen Baek, and Nando Fioretto. \"Constrained synthesis with projected diffusion models.\" Advances in Neural Information Processing Systems 37 (2024): 89307-89333."}, {"review_id": "Sh4YY3GLG2", "reviewer": "Reviewer_FfpW", "summary": "This paper presents Physics-Constrained Flow Matching, a method that enables arbitrary constraint enforcement on generative flow matching outputs through inference time constraint imposition. PCFM applies constraints by projecting the estimated final state $u_1$ (which is computed using a forward ODE solver) and then reversing the ODE to compute the next state $u_{\\tau'}$. As evaluated on a series of experiments evaluating performance on challenging PDE dynamics, PCFM outperforms baselines in qualitative metrics while ensuring perfect constraint adherence.", "questions": "Questions:\n1. Has any analysis been taken for out-of-distribution constraints (e.g., constraints that aren't necessarily satisfied by the training data)?\n2. Can the authors elaborate further on the SMSE performs worse as the number of flow matching steps increases? It seems intuitive that the MMSE decreases as the number of steps increases, but I'm having trouble reasoning as to its impact on SMSE.", "limitations": "The primary limitation that seems to be missing is an explanation of the impact of PCFM on runtime. It would be helpful if the authors could comment on this point.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths:\n1. **General Applicability (Significance):** The method provides several important selling points: *(1) Arbitrary Constraint Enforcement:* PCFM can handle arbitrarily complex constraint sets, a property which is demonstrated on complex PDE constraint enforcement, *(2) Inference-Time Method:* PCFM does not require retraining, so can be easily adopted by existing pretrained models.\n2. **Strong Empirical Performance (Significance):** PCFM reports strong results on the settings explored. The experimental settings are interesting, including varying degrees of constraint complexity on several different settings.\n3. **Writing and Presentation (Clarity):** The paper is well written and easy to follow. The method is presented clearly and well described.\n\nWeaknesses:\n1. **Inadequate Comparison to Similar Work (Originality, Quality):** The use of intermediate projections for continuous-time generative models is not a foreign concept in constrained generation literature [1-2]. For example, [2] applies projections to intermediate states in the diffusion reverse process, providing many of the same benefits that PCFM boasts (e.g., arbitrary constraint satisfaction, zero-shot constraint imposition). Furthermore, as a projection operator has already been designed for the PCFM's inference, it would be simple to add [2] as a baseline in the experiments. Due to the overlap between the methodology in these papers, a comparison would help clarify the novelty and performance claims.\n2. **Runtime Analysis (Significance):** PCFM relies on additional flow matching steps to balance constraint correction and sample quality. This (along with the complexity of computing the projection operator) likely contribute to much higher runtimes. Additional analysis should be included comparing the runtime of PCFM to the baselines. While the analysis in Appendix L.2 is interesting, it could be further expanded to include these computational overhead benchmarks as well.\n3. **Generality of PCFM (Significance):** The experimental analysis focuses on synthetic data for PDE governed dynamics. While synthetic data is generally sufficient for papers focusing on PDEs [3-4], PCFM is presented as a general method (e.g. lines 14-16, lines 67-71) for constraint satisfaction in other settings as well. These claims would be better supported by including \"real-world\" experiments to supplement the synthetic settings examined.\n\n\n---\n\n[1] Cheng, Chaoran, et al. \"Gradient-free generation for hard-constrained systems.\" (2025).\n\n[2] Christopher, Jacob K., Stephen Baek, and Nando Fioretto. \"Constrained synthesis with projected diffusion models.\" Advances in Neural Information Processing Systems 37 (2024): 89307-89333.\n\n[3] Huang, Jiahe, et al. \"DiffusionPDE: Generative PDE-solving under partial observation.\" arXiv preprint arXiv:2406.17763 (2024).\n\n[4] Ben-Hamu, Heli, et al. \"D-flow: Differentiating through flows for controlled generation.\" arXiv preprint arXiv:2402.14017 (2024).", "was_revised": true, "final_justification": "The authors have successfully addressed my primary concerns during the rebuttal period. To summarize:\n\n- **Comparison to ECI and PDM:** The additional results and explanations provided have addressed my primary concerns. Their inclusion in the final manuscript will mitigate this issue. While I believe there are some remaining questions of novelty (as the authors even state that their work is a generalized version of ECI), these additional points assist in differentiating the work from baselines.\n\n- **Runtime:** This should be acknowledged in the limitations. The table provided should be added to the paper. It's inclusion addresses this concern.\n\n- **Generality Claims:** The authors state their intent to clarify this and soften the claims in the updated manuscript. This will address my concern.\n\n- **Out-of-Distribution Settings:** I believe this is interesting, and its inclusion would further strengthen the analysis.\n\nOverall, while I still maintain some concerns about novelty, I have raised my score from a 3 to a 4 as I believe these clarifications push the paper over the acceptance threshold.", "rebuttal": "We thank the reviewer for the thoughtful assessment, noting the generality of our method and clarity of presentation. We appreciate recognition of PCFM’s handling of complex constraints at inference and strong empirical performance on challenging PDE benchmarks.\n\n## R3-1: Comparison to Prior Work on Intermediate Projection \n\nWhile PCFM shares inference-time constraint satisfaction goals, it differs from ECI and Projected Diffusion Models (PDM) in generality, theory, and extensibility.\n\n### Comparison to ECI:\n    \n1. ECI is a special case of PCFM (Appendix F), limited to linear, uncoupled constraints without correction penalty (λ = 0). PCFM generalizes this to support nonlinear and coupled constraints, such as nonlocal mass conservation.\n\n2. ECI is tested on simple linear cases and lacks scalability for nonlinear/global constraints. Its correction involves a full constrained optimization, which is expensive in large systems. PCFM uses lightweight projections during generation for guidance and final correction for feasibility.\n\n3. Though described as \"gradient-free\", ECI relies on analytical projections, limiting extensibility as robust optimization methods are gradient-based methods, which makes the extension to nonlinear constraints nontrivial. PCFM uses our custom batched differentiable solver for scalable PDE enforcement.\n\n### Comparison with PDM\n   \nWhile PDM also uses projection during sampling, its assumptions and modeling context differ fundamentally from PCFM.\n\n1. PDM enforces constraints at every timestep, which can over-constrain the trajectory, especially under nonlinear or global constraints like ICs or mass conservation that only apply at the final state. PCFM avoids this by guiding the trajectory through reversible integration and applying a single final correction, preserving consistency and flexibility.\n\n2. PDM is built on stochastic diffusion models, while PCFM uses deterministic flow matching. Though related via the probability flow ODE \\[10], their sampling behaviors and practical implementation differ.\n\n3. PDM has not been evaluated on PDEs with shocks or discontinuities, and its scalability to nonlinear global constraints is unclear. In contrast, PCFM is designed for such regimes and shows strong performance on complex PDEs like Burgers and Navier–Stokes.\n\nTo ensure fair comparison, we implemented a PDM-style ablation on our flow models. While it performs comparably on simpler cases like Heat, PCFM clearly outperforms on harder nonlinear dynamics, validating our design choices.\n\n### Table A: Comparison between PDM and PCFM\n| Dataset          | Method | MMSE      | SMSE      | CE (IC)      | CE (CL)  |\n|------------------|--------|-----------|-----------|--------------|----------|\n| Heat             | PDM    | 0.0045    | **0.0002**   | 0   | 0|\n| Heat    | PCFM      |    **0.00241**    |   0.00937     |         0    |       0       |\n| Burgers          | PDM    | 1.5309    | 0.0262    | 0  | 0|\n| Burgers | PCFM      |   **0.00052**        |     **0.00272**      |      0       |      0         |\n| Reaction-Diffusion | PDM    | 0.0174    | **0.0043**    | 0   | 0|\n| Reaction-Diffusion      | PCFM      |    **0.00026**       |     **0.00583**      |     0       |        0       |\n| Navier-Stokes    | PDM    | 0.1395    | 0.0786    | 0  | 0|\n| Navier-Stokes    | PCFM   | **0.0459**    | **0.0417**    | 0            | 0        |\n    \n\n## R3-2: Runtime Analysis and Overhead\n\nWe agree that runtime is important for assessing PCFM’s practicality. Accordingly, we include detailed runtime and memory comparisons using consistent batch sizes (8 for Navier–Stokes, 32 otherwise) and 100 flow steps on a single NVIDIA V100 GPU.\n\nPCFM incurs a modest overhead relative to unconstrained baselines like vanilla FFM, due to exact constraint enforcement. However, it achieves 4–6× speedups over DFlow \\[4], which requires costly backpropagation through the full ODE unroll. DiffusionPDE is faster but struggles with hard constraints, especially for nonlinear problems like Burgers and Reaction-Diffusion. ECI’s runtime grows with mixing steps but still fails to satisfy nonlinear constraints.\n\nFor linear problems like Heat and Navier–Stokes, projections converge quickly, reflected in PDM’s runtime. Overall, PCFM adds a small, necessary overhead for robust constraint enforcement. We will summarize these findings in the main text.\n\n### Table B: Runtime and Memory Comparison\n\n| **Dataset**          | **Metric**       | **PCFM** | **ECI** | **DiffusionPDE** | **D-Flow** |  **PDM**  | **Vanilla** |\n| -------------------- | ---------------- | :------: | :-----: | :--------------: | :--------: | :-------: | :---------: |\n| **Heat**    | Time (ms/sample) |  291.8  |  65.6  |      128.6      |   2770.5  |   361.3  |  65.2  |\n|                      | Peak Memory (GB) |   2.41   |   1.21  |       2.39       |    3.58    |    1.06   |   1.21  |\n| **Burgers** | Time (ms/sample) |  1371.0 |  780.6 |      211.2      |   8048.9  |   444.9  |  105.4 |\n|                      | Peak Memory (GB) |   2.57   |   1.26  |       2.43       |    3.65    |    1.08   |   1.22  |\n| **Reaction-Diffusion** | Time (ms/sample) | 636.6 | 744.4 | 159.9     |  6578.7 |   474.7  |  80.7 |\n|                      | Peak Memory (GB) |   2.97  |  1.72  |     2.93       |   4.33    |    1.43 |   1.43  |\n| **Navier–Stokes**    | Time (ms/sample) | 350.7 | 3355.6 | 678.8      | 16447.2   | 338.7  | 343.8 |\n|     | Peak Memory (GB) | 7.39  |   2.03    |  3.90     |  5.85  |  1.93   | 2.02 |  \n\n## R3-3: Generality Beyond Synthetic PDE Data\n\nWe appreciate the reviewer’s suggestion and will clarify the intended scope of our generality claim. Specifically, by \"general\", we refer to PCFM’s ability to support **a wide variety of constraint types**—including **nonlinear**, **coupled**, and **global equality constraints**—on a **general class of widely studied PDE systems**. We chose flow matching models (FFMs) to reflect the functional nature of PDE data [1].\n\nTo ensure comparability with prior work, we focused on **standard benchmark PDE datasets**, including **smooth equations** such as Heat and Reaction–Diffusion, and **nonlinear, shock-forming systems** such as Burgers and Navier–Stokes. These settings are commonly used in the constrained generative modeling literature \\[2, 3, 4], neural operator learning \\[5, 6], and foundation PDE models \\[7, 8].\n\nWe agree with the reviewer and will **soften claims about real-world generality** in the paper. While our method is broadly applicable to other constrained generation tasks (e.g., control, imaging, or molecular design), we reserve such extensions for future work.\n    \n## R3-4: Q1\n\nAs per the reviewer's request, we also include an out-of-distribution (OOD) constraint experiment to assess PCFM's generalization beyond constraints implicitly present in the data during training.\n    \nSpecifically, we consider **initial condition (IC) shifts** in the **Burgers equation**, where the pretrained model has not seen such constraints during training. We generated a new test dataset where the location of the initial sigmoid profile $u(x,0) = \\frac{1}{1 + \\exp\\left(\\frac{x - p_{\\text{loc}}}{\\epsilon}\\right)}$ is sampled from the intervals $[0.0, 0.1]$ and $[0.9, 1.0]$, while the training set used $p_{\\text{loc}} \\in [0.2, 0.8]$. The pretrained FFM model is **reused without retraining**, and all sampling methods are evaluated on the same OOD dataset to ensure fair comparison.\n    \nWe adopt 200 flow matching steps for the Burgers problem and the table below summarizes the results across methods.\n    \n### Table C: Performance on out-of-distribution dataset\n\n| Method       | MMSE       | SMSE       | CE (IC)     | CE (CL)     |\n| ------------ | ---------- | ---------- | ----------- | ----------- |\n| **PCFM**     | **0.0183** | **0.0140** | **0**       | **0**       |\n| ECI          | 0.1330     | 0.0565     | 0           | 3.6443e+00  |\n| DFlow        | 0.1559     | 0.0782     | 6.0027e+00  | 9.2038e-02  |\n| Vanilla      | 0.1884     | 0.0788     | 6.3685e+00  | 6.9078e-02  |\n| diffusionPDE | 0.1949     | 0.0803     | 6.4240e+00  | 6.2256e-02  |\n| PDM          | 2.0606     | 0.0281     | **0**       | 0           |\n\nPCFM achieves the lowest MMSE and SMSE while satisfying both ICs and nonlinear mass conservation even in this challenging OOD regime. These results demonstrate the **robust generalization** of our inference-time constraint enforcement, without requiring retraining or dataset-specific tuning. We will include this OOD evaluation in the final version of the paper.\n\n## R3-4: Q2\n\nWe thank the reviewer for the insightful question. As shown in Figure 9, increasing the number of flow matching steps improves MMSE by enabling the learned vector field to guide samples more accurately toward high-likelihood regions. However, we observe a mild increase in SMSE, particularly when stronger penalty weights are applied.\n\nThis can be attributed to reduced sample diversity: more frequent projection-based corrections during integration may cause samples to collapse toward a narrower region of the constraint manifold, lowering ensemble spread. This behavior reflects a trade-off between sample fidelity (MMSE) and diversity (SMSE).\n\nFlow steps and step size are user-controlled hyperparameters in flow matching frameworks [9]. In our case, the flow integration length may slightly suppress SMSE due to stronger adherence to the constraint manifold. However, this is tunable and can be adapted based on application needs.\n\n[1] Kerrigan et al., arXiv:2305.17209, 2023.\n\n[2] Cheng et al., 2025.\n\n[3] Huang et al., arXiv:2406.17763, 2024.\n\n[4] Ben-Hamu et al., arXiv:2402.14017, 2024.\n\n[5] Li et al., arXiv:2010.08895, 2020.\n\n[6] Li et al., arXiv:2111.03794, 2021.\n\n[7] Herde et al., arXiv:2405.19101, 2024.\n\n[8] Sun et al., arXiv:2404.12355, 2024.\n\n[9] Lipman et al., arXiv:2412.06264, 2024.\n\n[10] Albergo et al., arXiv:2303.08797, 2023."}, {"review_id": "7EedTAjW4a", "reviewer": "Reviewer_5iKF", "summary": "This paper proposes a method to guide the inference-time generation process of a flow-based vector field model $v_\\theta$, which has been trained to simulate the behavior of physical systems governed by partial differential equations (PDEs), so that the outputs satisfy physical constraints such as initial/boundary conditions and conservation laws (both linear and nonlinear), without requiring retraining of $v_\\theta$.\nTo generate the state $u_\\tau$ for $\\tau \\in [0,1]$, the method first solves an ODE based on $v_\\theta$ forward from $u_\\tau$ to estimate $u_1$. The estimated $u_1$ is then projected onto the manifold defined by the given constraints using a Gauss–Newton projection. An ODE is then solved backward from the projected $u_1$ to determine a correction target for $u_\\tau$. Finally, a penalized correction (Eq. (6)) based on this target is applied to obtain the final corrected $u_\\tau$.\n\nThe proposed method is applied to functional flow models (FFMs) trained to simulate physical systems such as the Navier–Stokes and Burgers’ equations. The results show that the method produces outputs that numerically satisfy the imposed constraints with high accuracy. In contrast to existing approaches, which often struggle to strictly enforce nonlinear constraints, the proposed method successfully generates solutions that comply with them. Furthermore, enforcing these constraints improves performance on evaluation metrics such as mean squared error (MSE) and FPD.", "questions": "### Comments and Typos\nIt seems that an additional operation such as $u \\leftarrow u_{\\tau'}$ may be required after line 8 in Algorithm 1.\n\nAlgorithm 1 should be explicitly referenced and discussed in the main text.\n\nLine 127: constrained -> Constrained?\n\nLine 136: $\\mathcal{F}(u; \\phi)$ -> $\\mathcal{F}_{\\phi}(u)$?", "limitations": "See weaknesses above.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 4, "originality": 2, "strengths_and_weaknesses": "### Strengths\n* The proposed method successfully imposes both linear and nonlinear physical constraints on flow-based generative models of physical systems in a zero-shot manner.  \n* The experiments demonstrate that the proposed method can generate solutions that satisfy the constraints, and that enforcing such constraints also improves the accuracy of the generated results.\n\n### Weaknesses\n* The proposed method shares a broadly similar procedure with the existing ECI framework, which limits its technical novelty. As far as I understand, the primary distinction lies in the correction step introduced in Eq. (6).\n* The experiments are limited to relatively small-scale settings with only two degrees of freedom.", "was_revised": true, "final_justification": "The authors' response clarified several aspects of the proposed method. While the technical novelty may be somewhat limited given its similarity to ECI, I find that the paper still makes a meaningful contribution. In particular, the experimental results demonstrate that PCFM handles nonlinear constraints more efficiently than existing methods, which I consider a valuable empirical finding.\n\nI have no major remaining concerns, and I maintain my original score in support of acceptance.", "rebuttal": "We thank the reviewer for their thoughtful and constructive review. We appreciate your recognition of the method’s ability to enforce both linear and nonlinear physical constraints on flow-based generative models in a **zero-shot** manner and your general positive inclination towards our method. Below, we address the specific concerns raised.\n\n## R4-1: Novelty Relative to the ECI Framework\n\nWe appreciate the reviewer’s observation regarding the similarities to the ECI framework. While our method shares a conceptual goal: enforcing constraints at inference time, it differs significantly in both **generality** and **mechanism**, and we clarify these distinctions below:\n\n* **Generality of Constraints**: ECI introduces a novel mixing-based correction strategy but is primarily demonstrated on **simple, linear, and decoupled constraints**. It does not provide a general framework for handling **nonlinear, coupled, or global constraints**, nor does it show robust performance on PDEs with **shocks or discontinuities**. In contrast, PCFM is **built to support arbitrary nonlinear equality constraints**, including **nonlocal mass conservation**, and shows strong empirical performance in such regimes.\n\n* **Theoretical Connection**: As detailed in Appendix F, ECI is a **special case** of PCFM when constraints are linear, non-overlapping, and the penalty is omitted (λ = 0). PCFM generalizes this setting to support more **complex and coupled constraints**.\n\n* **Extensibility and Differentiability**: Although ECI is described as “gradient-free,” it depends on **analytical projections** that do not extend naturally to nonlinear or global constraints. PCFM uses **general-purpose automatic differentiation** and **Jacobian-vector products**, enabling scalable enforcement in **high-dimensional and nonlinear settings**.\n\n* **Inference Strategy**: ECI corrects the final output through constrained optimization, which can be expensive in large systems. PCFM instead uses **efficient projection steps during generation** and applies a **single correction step** at the end (Eq. 6) to ensure feasibility without high overhead.\n\nWe agree that Eq. 6 is a key contribution: it ensures feasibility while promoting alignment with the learned flow field, preserving generative consistency. We validate this in Appendix L.2, where PCFM outperforms ECI even under linear constraints, and more strongly under nonlinear and coupled constraints.\n\nThat said, Eq. 6 is part of a broader pipeline. PCFM integrates projection-based solvers, reversible integration, and constraint-guided updates, drawing from numerical PDE solvers to enable general, zero-shot constraint enforcement. Unlike ECI, which solves a constrained optimization problem at every sampling step, PCFM uses lightweight intermediate projections and a final correction. To ensure scalability and efficiency, we implement a batched, differentiable projection solver using reverse-mode Automatic Differentiation and Schur-complement updates, avoiding retraining while supporting complex constraints.\n\n## R4-2: Experimental Scope — Small-Scale Systems\n\nWe respectfully clarify that while the reviewer refers to our setting as having “only two degrees of freedom,” the Navier–Stokes system is inherently three-dimensional, spanning two spatial dimensions and one temporal dimension. Moreover, we enforce three types of constraints simultaneously—initial conditions, local consistency, and nonlinear global conservation laws—highlighting the method’s applicability to high-dimensional and tightly constrained PDE systems.\n\nTo ensure thorough evaluation, we included PDEs of varying difficulty: from relatively simple systems like the Heat equation, to intermediate ones such as Reaction–Diffusion, and finally to nonlinear, shock-prone systems like Inviscid Burgers and Navier–Stokes.\n\nAll datasets and spatiotemporal resolutions are sourced directly from prior work (Appendix B). Evaluating on 2D and 3D PDE systems is standard practice in the literature on neural operator learning [1, 2], scientific foundation models [3, 4], and constrained generative modeling [5, 6], making our experimental setup both representative and widely accepted.\n\n\n[1] Li, Zongyi, et al. \"Fourier Neural Operator for Parametric Partial Differential Equations.\" NeurIPS (2020).\n\n[2] Li, Zongyi, et al. \"Physics-Informed Neural Operator for Learning Partial Differential Equations.\" arXiv:2111.03794 (2021).\n\n[3] Herde, Maximilian, et al. \"Poseidon: Efficient Foundation Models for PDEs.\" arXiv:2405.19101 (2024).\n\n[4] Sun, Jingmin, et al. \"Towards a Foundation Model for Partial Differential Equation: Multi-Operator Learning and Extrapolation.\" arXiv:2404.12355 (2024).\n\n[5] Huang, Jiahe, et al. \"DiffusionPDE: Generative PDE-Solving under Partial Observation.\" arXiv:2406.17763 (2024).\n\n[6] Négiar, Geoffrey, Michael W. Mahoney, and Aditi S. Krishnapriyan. \"Learning Differentiable Solvers for Systems with Hard Constraints.\" arXiv:2207.08675 (2022).\n\nThank you for pointing out other minor presentation issues. In the revised manuscript, we will fix those typographical errors."}]

Author rebuttal:

Meta-review: The paper presents a method termed Physics-Constrained Flow Matching (PCFM), which enforces physical constraints in flow-based generative models via inference time guidance. This is implemented using the Gauss-Newton method.

The reviewers agreed that the paper is well written and that it contains significant results: the method is generally applicable and demonstrates strong empirical performance.

There were some concerns about the term "hard constraints" and about the novelty compared to ECI abd PDM but they were addressed in the rebuttal.

After the rebuttal, all reviewers recommend acceptance, provided that the authors update the paper according to the discussion (including relaxing the term "hard constraints").
Decision: Accept (poster)