Learning Reduced Fluid Dynamics
Authors: zherong pan
Xifeng Gao
Kui Wu
Venue: Submitted to ICLR 2023 (2023)
Keywords: Fluid Dynamics
Model Reduction

Abstract: Predicting the state evolution of ultra high-dimensional, time-reversible fluid dynamic system is a crucial but computationally expensive task. Model-reduction has been proven to be an effective method to reduce the computational cost by learning a low-dimensional state embedding. However, existing reduced models are irrespective of either the time reversible property or the nonlinear dynamics, leading to sub-optimal performance. We propose a model-based approach to identify locally optimal, model-reduced, time reversible, nonlinear fluid dynamic systems. Our main idea is to use stochastic Riemann optimization to obtain a high-quality a reduced fluid model by minimizing the expected trajectory-wise model-reduction error over a given distribution of initial conditions. To this end, our method formulates the reduced fluid dynamics as an invertible state transfer function parameterized by the reduced subspace. We further show that the reduced trajectories are differentiable with respect to the subspace bases over the entire Grassmannian manifold, under proper choices of timestep sizes and numerical integrators. Finally, we propose a loss function measuring the trajectory-wise discrepancy between the original and reduced models. By tensor precomputation, we show that gradient information of such loss functions can be evaluated efficiently over a long trajectory without time-integrating the high-dimensional dynamic system. Through evaluations on a row of simulation benchmarks, we show that our method lower the discrepancy by 45%-97% over conventional reduced models.

Under review as a conference paper at ICLR 2023
LEARNING REDUCED FLUID DYNAMICS
Anonymous authors
Paper under double-blind review
ABSTRACT
Predicting the state evolution of ultra high-dimensional, time-reversible fluid dy- namic systems is crucial but computationally expensive task. Model reduction has been proven an effective method to reduce computational costs by learning low-dimensional state embedding. However, existing reduced models are irre- spective of either the time reversible property or the nonlinear dynamics, leading sub-optimal performance. We propose a model-basee approach to identify lo- cally optimal, model-reduced, time reversible, nonlinear fluid dynamic systems. Our main idea to use stochastic Riemann optimization to obtain high-quality reduced fluid model by minimizing the expected trajectory-wise model reduction error over a given distribution of initial conditions. To this end, our method for- mulates the reduced fluid dynamics eterized by the reduced subspace. We further show that the reduced trajectories are differentiable with respect to the subspace bases over the entire Grassman- nian manifold, under proper choices of timestep sizes and numerical integrators. Finally, we propose loss function measuring the trajectory-wise discrepancy between the original and reduced models. By tensor precomputation, we show that gradient information of such loss function can be evaluated efficiently over a long trajectory without time- -nteinnggga the high-dimensional dynamic system. Through evaluations on row of simulation benchmarks, we show that our method reduces the discrepancy by 50% - 90% over conventional reduced models.
invertible state transfer function param-
INTRODUCTION
High-dimensional Partial Differential Equations (PDE), especially fluid dynamic systems, find vast applications in the field of scientific computation Moin & Mahesh (1998); Alfonsi (2009), PDE- constrained optimization Biegler et al. (2003); Herzog & Kunisch (2010), design prototyping Baysal & Eleshaky (1992); Zang & Green (1999), fluidic devices design Du al. (2020); Li et (2022), and digital entertainment Bridson & Batty (2010); Bridson (2015), to name a few. A fundamen- tal task all these applications lies in the efficient prediction of numerical solutions over a long horizon. In design prototyping, for example, designer needs to quickly preview the fluid flow surrounding an acrial vehicle in order to refine its form factor. In a game engine, a fluid simulator needs to achieve real-time performance provide interactive special effects for players. Although abundant numerical tools Petrila & Trif (2004); Demkowicz et al. (1989) have been developed over the past decades with improved efficacy, their algorithmic complexity still challenging the limits of current computational resources. As parallel effort, the idealized, incompressible, inviscid Eu- lerian fluid should be time reversible and energy preserving Duponcheel al. (2008), and dedicated numerical schemes are proposed to faithfully preserve these properties in discrete setting Row- ley & Marsden (2002); Pavlov al. (2011). This implies that the initial condition of trajectory can be recovered from any state thereafter and the discrete total energy is constant throughout the predicted trajectory. Although idealized fluid models are not pursued applications, their accurate prediction is an important criterion of reliable numerical schemes.
Since their proposal Berkooz et al. (1993); Rowley (2005), model reduction has been quickly estab- lished as one of the most effective approaches that can significantly reduce the PDE prediction cost. By restricting the state variables to dinedsmiiionnl linear and nonlinear mani-olssssds, the dimen- sion of associated dynamic system be reduced by orders of magnitude. Over the years, several data-driven and data-free approaches have been proposed to identify sub-manifolds that can capture the complex dynamic behaviors of fluids. The earliest data-driven method of Proper Orthogonal

---

Under review as a conference paper at ICLR 2023
Decomposition (POD) Berkooz et al. (1993) finds the optimal linear subspace that best explains the variation of the state distribution. However, POD is flawed in that it ignores the temporal dependence of state variables. This problem is remedied by the Dynamic Model Decomposition (DMD) Schmid (2010) that finds the optimal linear subspace that best approximates the Koopman operator. How- ever, these data-driven algorithms are irrespective of the nonlinearity in the underlying PDE. Com- paratively, data-free methods, such as balanced POD Rowley (2005), Hl--optimization Gugercin al. (2006), and modal analysis Taira al. (2017), identify bases corresponding the intrinsic property of PDE by analyzing the system transfer matrices in the frequency domain, and are thus independent of data. Unfortunately, these techniques are largely limited linear systems and their extensions to nonlinear fluid dynamics, such as Yang et al. (2019), are in their infancy.
More generally, the construction of reduced fluid models has been formulated as machine learning problems for system identification. The vast majority of prior works generalize the non-intrusive approach and identify the state transfer function via supervised learning in an existing sub-manin man-manifold, where the transfer functions are parameterized via radial basis functions Zhang et al. (2016), feed- forward networks Hsich Tang (1998), recurrent networks Pearlmutter (1989); Wang et al. (2018), etc. More recent approaches jointly learn the state transfer function and identify the sub-manifold via convolutional autoencoder Wu et al. (2021); Hasegawa et al. (2020). Unfortunately, all these non-intrusive learning techniques cannot preserve the time reversible property of idealized fluid, potentially leading to large prediction error or requiring a large dataset.
Figure 1: Given a distribution of initial conditions I, we identify reduced-order fluid model 5' (e, U) by optimizing the bases that minimize the expected trajectory -wite discrepancy loss Cayn- Our output model e² (v, U) can perform efficient and as as-accurate- as-possible fluid trajectory predictions.
8t (TP .))
(66600)
(""" .))
Caya
Caya
We propose a machine learning approach to identify locally optimal, time reversible, reduced- -odder fluid dynamic models. We first interpret the linear subspace of fluid velocities as point on the Grassmannian manifold and study the dependence of reduced trajectories on the choice of subspace. Thanks to the time reversibility, we show that the map from the subspace bases to reduced trajec- tories is globally differentiable, which allows us to optimize the reduced model via gradient-based Riemannian optimization. We further propose a trajectory-wisc discrepancy loss that penalizes the difference between the full-order and the reduced trajectories. To make the optimization tractable, we propose tensor precomputation scheme to accelerate the pooo-propagttion of gradient informa- tion. Figure 1 illustrates the high-level pipeline of our method that fine-tunes the reduced fluid model to minimize the expected trajectory-wee wise discrepancy loss over the distribution of initial conditions. In essence, our method extends prior optimal reduced bases construction algorithm Berkooz et al. (1993); Schmid (2010) to the nonlinear, idcalized fluid dynamic model. As an intrusive approach, method preserves the desirable property of time reversibility. When compared with POD-type reduced model baseline on a row of idealized fluid simulation benchmarks, our method lowers the discrepancy by 50% 90%.
2 RELATED WORK
We review related works machine learning for solving ODE and PDE, reduced physics models beyond fluid dynamics, and finally learning under hard constraints.
Learning for Solving ODE and PDE:
To study the complex behavior of dynamic systems, var- ious mathematical models have been proposed for idcalized models of fluid, solid, elasto-magnetic fields, etc. However, there are oftentimes subtle discrepancies between these models and real-world observations that are hard to model, in which cases machine learning stands out as an effective ap- proach for acquiring these behaviors from groundtruth data. Chen et al. (2018) propose to learn such dynamics as a general Ordinary Differential Equation (ODE) with the time derivative of state predicted via a neural network. Although this method is applicable to general dynamic systems, it does not reflect the spatial and temporal structures of certain systems, which limits its accuracy, data-efficacy. and scalability to sigh--imensional systems such as fluids. Several follow-up works improve the network architecture to reflect additional structures. For example, the inter-dependency
2

---

Under review as a conference paper at ICLR 2023
among spatial variables is oftentimes local and sparse, which could be modeled via neighborhood message passing Battaglia et al. (2016); Li et al. (2019). Hamiltonian dynamics are time reversible and energy preserving, which is modeled by learning the Hamiltonian operator in "canonical" co- ordinates Greydanus et al. (2019), generalized coordinates Cranmer et al. (2020), or ambient space with additional constraints Finzi et al. (2020). However, the above techniques are using La- grangian coordinates, while fluid mechanics are oftentimes modeled via an Eulerian grid, see e.g. Takahashi et al. (2021), which is a major point of difference from our method. Parallel efforts have been made to learn Eulerian fluid mechanics Um et al. (2020); Takahashi et al. (2021); Holl et al. (2020); Prantl et al. (2019); Kim et al. (2019). Some of these works Um et al. (2020); Takahashi et al. (2021); Holl et al. (2020) learn to control fluids via differentiable sim- ulators but the dynamic systems are not learned. Other works Prantl et al. (2019); Kim et al. (2019) learn to predict short future trajectories of free-surface flows. As the major difference from these techniques, our goal is to predict arbitrarily long trajectories by utilizing the time reversible structure of the dynamic system to guarantee stability. On the downside, however, our method cannot predict free-surface flows
Model reduction is a special kind of dimension reduction technique dealing with time series datasets and we refer readers to Rowley & Dawson (2017) for a review of its application in fluid mechanics. Other than fluid, reduced models have been adopted in predicting the behaviors of solid Sampaio & Soize (2007), electromagnetic fields Ralph-Uwe et al. (2008), quantum and molecular mechanics Mohan & Fredrickson (2020), neuron propagations Am- sallem & Nordstrom (2016), etc. A successful reduced model involves two steps: 1) embed the data into a proper subspace that well explains the data variations; 2) project the dynamic system into the subspace. Conventional techniques for model reduction are restricted to linear dynamic systems, for which optimal linear subspace can be identified via POD or DMD Berkooz et al. (1993); Rowley (2005) and the projected dynamic system can be precomputed via Galerkin projection. More gen- eral machine learning techniques have been proposed for an extension to nonlinear dynamics. For example, convolution autoencoder has been used to identify nonlinear subspaces Wu et al. (2021); Hasegawa et al. (2020). The ROM-net Daniel et al. (2020) learns to select a suitable subspace from a dictionary. Li et al. (2017) proposes to represent the linear subspace bases the output of univer- sal neural network. In order to efficiently project the nonlinear dynamic system into the subspace, the Discrete Empirical Interpolation Method (DEIM) Chaturantabut & Sorensen (2010) has been proposed to select a sparse set of interpolation points. The interpolation points are then contracted with the subspace bases in an intrusive manner. Non-intrusive approaches use universal neural net- works to learn the entire nonlinear transfer function Wu et al. (2021); Hasegawa et al. (2020); Lee al. (2021) or part of the nonlinear terms Maulik et al. (2019). It has been noticed in Amsallem & Nordstrom (2016) 6); Liu et al. (2015) that time reversibility and energy preservation features can be preserved by using an intrusive approach, which a main reason behind our technical choice.
Learning Reduced Physical Models:
Learning Under Hard Constraints: Our work deals with idealized fluid satisfying two hard con- straints: 1) incompressibility and 2) time reversibility. Since prominent training algorithms Duchi et al. (2011); Kingma & Ba (2014) and neural network architectures are designed for unconstrained optimization, dealing with hard constraints has been a long-standing problem Márquez-Neila et al. (2017). There are two general approaches to inform a learned model of hard constraints: soften- ing and constraint layers. Softening transforms the hard constraint into soft losses and relies on unconstrained optimizations. Some hard constraints model invariant variables, in which case data augmentation could be used to enforce a neural network gives the same output over all invariant transforms of inputs. In the learning of physical models, softening has been adopted enforce physical correctness Sirignano & Spiliopoulos (2018); Ober-Blöbaum & Offen (2022), fluid incom- pressible Ajuria Illaramendi et al. (2020), and collision-fre free constraints Tan et al. (2022), and data augmentation has been used to enforce invariance to rigid Morozov et al. (2021) and Galilean trans- formations Ling et al. (2016). A common problem with all these approaches lies in the unpredictable constraint violation in regions of insufficient data coverage. To exactly impose hard constraints, a series of works Amos & Kolter (2017); Agrawal et al. (2019) propose to formulate the constrained optimization as differentiable layer in the neural network architecture. In particular, the entire fluid simulator has been formulated as a differentiable layer Schenck & Fox (2018); Takahashi et al. (2021) for model- based control and system identification. The incompressible constraint has also been formulated as an elliptic PDE solver layer in Mohan et al. (2020). Although these techniques can enforce hard constraints, the cost of forward- and hack-propagatioos through these layers
3

---

Under review as a conference paper at ICLR 2023
prohibitive. Even worse, these layers must be evaluated during both training and test time. Our method uses the constraint layer approach to enforce fluid incompressibility and time reversibility, by incorporating the reduced model Liu et al. (2015) our differentiable layer. However, code the constraint property into the reduced bases, which is fixed during test time, leading to the low computational cost of trajectory prediction.
TIME REVERSIBLE REDUCED FLUID MODEL
We briefly review the underlying geometric structure and associated computational model of ideal- ized, incompressible, inviscid fluid Pavlov et al. (2011). Given a simulation domain M, the fluid configuration can be described as vector field v ((M) where v(z) for any x . represents the velocity of fluid at z. The governing equation for is:
V.
(I)
where À is the pressure field, which is also the Lagrangian multiplier for the divergence- preegencee-fre con- O. The above system is closed with appropriate initial and boundary conditions. Pavlov et al. (2011) proposed time-reversible, energy preserving spatial and temporal discretization schemes for Equation large-scale nonlinear system equations. Reduced-order model Liu et al. (2015) scales down the cost by embedding v into a p-dimensional subspace with divergence- free, orthogonal basis U, giving Us where is the coefficient vector. The reduced-order governing equation can be derived via Galerkin projection:
straint V
However, directly time integrating the discrete system requires solving
2+ Jm V (U2) =U =,
(2)
where the second term is the reduced-order advector, which could be succinctly written as a con- traction with third-order tensor Ckij:
2k ΣΣ Ckijziz = s.t. Cki SNL ( (Uk.VXU;xU)) i
where we use 21 (resp. Uk) to denote the kth element (resp. column). For fast reduced trajectory prediction, the tensor Ckij is precomputed, and small p is used. An essential feature of Cust is antisymmetry: Ckis -Cjik, which implies that the continuous-time, reduced system is also mmery--presenve as:
d
|2| =2 Eas E(Ckis CCak) 0.
dt
kij
Using a variational integrator, e.g. the trapezoidal rule, the energy will also be conserved in time-discrete computational model. We use a superscript to denote variables at the next time
instance, the superscript denotes the variable at the dth timestep, and òt denotes timestep size. The trapezoidal rule relates 2t and 2 by:
= + -(=))=0 x* x=))
s.t. =(+)) = C 2 +
2 2j
(4)
51
2
2
from which 2 can be solved via the Newton-Raphson Raphson method to satisfy 1|+ |2|2, i.e. energy conservation, as well as discrete-time reversibility. These remarked features, originally discovered in Pavlov et al. (201 I); Liu et al. (2015), achieve ideal balance between computational efficacy and numerical stability. As pointed out by Pavlov et al. (2011), although real-world flows are simulating ideal flows is crucial benchmark for evaluating the stability and fidelity of simulator. More general non-reversible flows can be modeled by adding additional constitutive terms. an example, we could add a viscous term pV . v) to model energy dissipation and this term can be projected to the reduced space via Galerkin projection. Figure 2, we plot the procedure both our learned reduced model and the groundtruth fullspace model. We formalize and prove these properties in Appendix A and Appendix and sufficiently small öt, we define the function 2* (2) by slight abuse of notations. The accuracy of reduced model relies on proper choice the basis vector U, which remains a difficult but underappreciated problem.
not ideally energy-presers
energy dissipation under different μ using
In particular, Equation 4 defines a unique 2* given
4

---

Under review as a conference paper at ICLR 2023
12.5
12
10.0
diffuse=0
10
diffuse=0
7.5
diffuse=1
diffuse=1
diffuse=10
diffuse=10
5.0
diffuse=100
diffuse=100
2.5
(a)
4
(b)
0.0
Figure 2:
We plot the energy dissipation cause by viscous term under 10, 100, simulated using our learned reduced model (a) and the groundtruth fullspace model (b).
4 REDUCED MODEL OPTIMIZATION
As illustrated in Figure 1, we propose to identify reduce-order fluid models via gradient-basee based op- timization of U to minimize the trajectory-wisst discrepancy between reduced-order model (Fqua- tion 3) and the full-order model (Equation 1). In this section, we first discretize the spatial computa- tional domain (Section 4.1), we then propose our discrepancy loss function (Section 4.3), and finally discuss optimization algorithm (Section 4.4).
4.1 SPATIAL DISCRETIZATION
We assume that M is discretized using tetrahedron mesh or rectangular grid via Discrete Exterior Calculus (DEC) as Pavlov et al. (2011). As a result, each vector field has a finite dimension p. We use a bar to denote discrete variable so ú R". U belongs the intersection of Stiefel manifold St p) and the divergence- free basis subspace: D(n.p) (Ü R"'P/V - = 0}, where V. (===m)xx is the discrete divergence operator and m >> is the dimension of divergence-free velocity subspace. The elements of U can also be identified with the elements of St(m, m). Indeed, we can find a set of unit, orthogonal bases D e IR nxm spanning the subspace of divergence- velocity fields. For each U, we identify some ÚTM St(m,p) p) such that U DUT As illustrated in Figure 8, a point on St(n,p) p) is the bases of a p-dimensional velocity field subspace, while a point on St m. p) the bases of p-dimensional divergence use U to project the velocity field into subspace, we are only interested in the Grassmannian Manifold (the manifold of velocity subspace irrespective of the particular bases), but we use Stiefel representation for better memory and computational efficacy. In other words, we treat U as our decision variable. We further write the tensor coefficient Ckt as a function C(U,, U,, U,), which is derived by discretizing the continuous definition of Ckij in Equation 3 using DEC.
velocity field subspace. Since we merely
4.2 LIFTING TRANSFER FUNCTION FROM REDUCED- TO FULL-SPACE
In order to optimize the accuracy of reduced dynamic system, we first need to compare simulated trajectories generated by different bases U. However, the coordinate vector 2 of different U is incomparable, as they reside in different linear subspaces. To resolve this problem, we propose to lift z to = Üz in the ambient space R", so that two vectors can be compared by the induced metric in the Euclidean space. Further, we can smoothly extend the reduced-order simulator function to the ambient space using the projection operator P U0 and P. = I:
J' = (0T
(5)
In other words, the velocity component orthogonal to the subspace is stationary, and the tangential velocity governed by the reduced dynamic system. As detailed in Appendix C, the above extension can be written as a function defined on the Grassmannian manifold: o (7, pT)) R2 x Gr(m. -
R", where we denote pm = - [ū"]". With the smooth extension, we can evaluate the derivatives of 5* with respect to E and the subspace. We can also compare two velocity fields generated by reduced-order simulators using different subspaces. Note the full-order dynamics (Equation 1) can be identified with " mmxm The above lifting is not unique, and a useful alternative is to discard the orthogonal component, i.e. setting P,"* which is discussed in Appendix C.2. As our major
5

---

Under review as a conference paper at ICLR 2023
contribution, show in Appendix C that the above function गt is well-defined smooth function on Gr(m,p), We further show that for any differentiable loss function L(v*), ). its derivatives with respect to the bases can be efficiently computed under proper representation of U as a manifold point.
Algorithm 1 Forward-Backwardd (o),)) (0)
Precompute tensor Ckij CÜk, U,) 2: Precompute tensor C(DD", ,0,0000000 U,)
3: for d= ".,,, T do
> forward propagation
4: jd+1 8 (vd (
5: G
backward propagation
6: Evaluate VL + 07 Lo"
35
7: for d = T 1,,,,,, do
8: G - G+ Equation 12
VL + ao
VL By (an
9:
O6d
10: Compute ViL via Equation 11
D divergence-free projection
Return vgL
Algorithm 2 RAMSGRAD(I,U)
Input: 32,32,
1: m + 0, 0, v 0, 0
2: while Not converge do
3: Sample n ~
>> we always use batch size equals to
4: g Barward-Backward(o" .0)
5:
m - --1 (1- B1)g
6:
61-+(1- (1-31))1 32)
7:
ú max(i,v)
8:
U Retract(U. (U, -amm
> by QR factorization approximate parallel transport
9:
P.m
Return Ú
4.3 REDUCED DISCREPANCY LOSS
The differentiable structure of reduced fluid allows us minimize the discrepancy between reduced- and full-order model in an efficient model-basse manner. Given two velocity fields and ot a full- order model should satisfy the governing equation of motion, which inspires the following discrep- ancy measure:
Laya (T*, ;0)
σ+ - D öt
DDT
C(DDT 5 +0 0- + +
66)
2
This is similar the physics correctness loss used in Sirignano & Spiliopoulos (2018); Ober- Blöbaum & Offen (2022) and we absorb the linear divergence-free constraint by using the projection operator DD Again, evaluating involves sparse linear solve for each of the 7 timesteps. But we can accelerate this computation thanks to the low-rank property of the velocity fields. Since, and 5t both reside in low-rank spaces, we can write:
T+ +
5 + ) C D,,Ū,)
2 + +
C(ŌDT
2
2
and precompute the tensor C(DDT U,, U,) via p² sparse linear solves at the cost of On For trajectory with T >> p² timesteps, this operator reduces the cost of evaluating Ldyn from (n")T) O(ni p Tnp?).

---

Under review as a conference paper at ICLR 2023
4.4 STOCHASTIC RIEMANN OPTIMIZATION
Using a low-dimensional subspace, it is impossible to approximate all fluid simulation trajectories with sufficient accuracy. Instead, reduced models are designed optimize subset of trajectories with a given distribution I of initial conditions, i.e. T" ~ I and our goal is to solve the following problem via stochastic Riemann optimization:
argmin Ego.I Σ Laya )
(eD(r.)))Ss(n..)))
where T is the horizon of trajectory and y E (0,1] is a constant discount factor. Riemann op- timization is a well-studied problem in both deterministic and stochastic settings and we use the RAMSGRAD algorithm proposed in Becigneul & Ganea (2019). This algorithm requires both the retraction and parallel transport operators on St(n.p). We use QR-factorization for the retraction operator Bendokat et al. (2020). Unfortunately, there is no efficient way to compute the parallel transport operator Edelman et al. (1998), so we approximate the transport operator by projecting out the non-tangential component. This corresponds to using a single step of forward Euler integrator to solve the associated ODE of the transport operator. Again due to time reversibility, the objec- tive function globally differentiable with respect to U under compact and sufficiently small öt. We outline our forward- seel-bacwwrtt gradient propagation in Algorithm 1 and adapted RAMSGRAD Algorithm 2. These algorithms are well-defined due to the following lemma:
Lemma 4.1. For any compact initial distribution I, there exists a sufficiently small öt, such that the objective function Σd-o 7at Ldm((T)) is globally differentiable, i.e. for any 2 I and Ü D(n.p) p) n pt..
Proof. Since I is compact, e is uniformly upper bounded by some and | ||O S r. By Corollary A.5, there exists sufficiently small öt making any 2d differentiable, reversible function of 2 This also implies Ed is a differentiable, reversible function of 7o under the definition of Equation 5, and our result follows.
EVALUATION
We implement our method using Pytorch with fluid simulator im- plemented via native C++ with CPU parallelization, and perform all the computations on an AMD Threadripper 3970X CPU having 32 cores. We initialize our method using a conventional POD-type algorithm. Given I, we first sample a set of N trajectories using the full-order dynamics (Equation 1) and then perform POD-type basis extraction. The number of extracted bases is determined by truncating the eigenvalues below of the largest eigenvalue. We always use a batch size of 1. The performance of our method is summarized in Table 1. We consider two variants of our method:
Time/ma)
2
10
40
50
Figure The cost of evaluating x x=()) plotted against p.
coupled case, where Caiz is treated as function C(U,, Uy) as discussed in Section 4, and de- coupled case, where Ckij is treated as an antisymmetric independent decision variable. Our main experiments are performed with the coupled case. Experiments with the decoupled case and sum- mary of decision variables are included in Appendix D. The efficacy of trajectory prediction using a reduced-orde -order model depends on p as illustrated in Figure 3, so the runtime performance of both the POD baseline and our method are the same, while the cost of evaluating the full-order model 252ms (26x slower than the reduced-order model with 49).
Benchmark
é- 0.05
c-
C- 0.001
< 0.0001
p Loss-POD Loss-Ours
Loss-POD POD
Loss-Ours
Loss-POI)
Loss: POD) Loss-Ours
p
Loss-POI> Loss-Ours
Taylor Vortices
8
4.84
0.58
11
3.93
0.39
16
1.99
0.17
25
0.70
0.07
Plume Rise
57.04
6.37
28.23
5.38
18.53
2.30
26
5.96
1.44
Plume Rise+Obst
133.16
10.48
46.47
8.60
16
17.57
2.57
30
5,71
1.05
Spherical Plume
36
120.89
44.89
Two Plume
59
103.23
49.22
Table 1: Summary of benchmarks for comparing POD and our method under different and
7

---

Under review as a conference paper at ICLR 2023
POD
coupled(8)
Our Method
coupled(11)
Discrepancy
Discrepancy
coupled(16)
coupled(25)
10
155
20
25
10 Elapsed Time (hr) (c)
20
#Bses
(b)
Figure 4: (a) Velocity magnitude field snapshots the Taylor vortices benchmark, generated by full-order model (top row), our method with e 0.0001 and p 25 (middle row), and POD with e p= 25 (bottom row). (b) Trajectory-wise discrepancy loss of POD and our method, under different p. (c) The convergence history of our method over 24 hours.
0.0001 and
Our first benchmark is Taylor vortices Pavlov et al. (2011), where two vortices are separated by a distance slightly larger than the critical threshold. We use velocity field discretized a 64 64 rectangular grid with the periodic boundary condition, leading to n = 8192. This is single trajectory (I is deterministic) and we set T 500. 86 0.01. We experiment with four parameters 0.01, 0.001, and 0.0001 and the number of bases is p 8, 11, 16, and 25, correspondingly. With each U as the initial guess, wise discrepancy loss against the number of bases p and the convergence history of our method. Compared with POD bases, our method reduces the discrepancy loss by 87.93%, 90. 90.12%. 91.47%, and 90.16%, respectively. Snapshots are shown in Figure 4a, where our method predicts velocity field closer to the full-order groundtruth.
0.05,
run our optimizer for 24 hours. In Figure 4bc, we plot the trajectory-
Our second benchmark involves having a smoke plume rise at a constant speed. We use rectangular domain of [0, 1 with all Dirichlet boundary conditions. The region of [0.25, 0.75] [0.125, 0.375] is occupied by the smoke with constant speed (0, 1) the remaining regions have zero velocity, and we use T = 1000. All other settings are the same as our first benchmark. The discrepancy loss and convergence history are plotted in Figure 13bc. We experiment with four parameters = 0.05, 0.01, 0.001, and 0.0001, the corresponding numbers of bases p are 6, 9, 15, and 26, respectively. Our method reduces the discrepancy loss by 88. 18.8 80.94%, 87.60%, and 75. 79%, respectively. We have also tested a variant of our method with an obstacle in the simulation domain, where our method reduces the discrepancy loss by 92. 81.49%, 85.38%, and 81.70% respectively. Snapshots of our second benchmark are shown in Figure 13a and Figure 14 of Appendix H.
20
coupled
(oupled(8)0.000
Loss-Init.
Loss-8k-Iter
Loss-12k-Iee 12k-Iter.
1.5
0.01
4.84
0.79
N/A
coupled(
(NO.5) couple di noise)
0.05
5.04
0.83
N/A
0.25
10.16
1.23
N/A
0.5
22.77
1.51
1.08
2500
5000
7500 llteration
10000 12500
Noise
21.06
2.37
N/A
Figure 5: The convergence history of four instances of learning reduced Taylor vortices with e = 0.05, p= 8, and difference noise levels ē= 0.01, 0.05, 0.25, 0.5. We first run the four training instances for 8000 iterations, which alrcady brings the ultimate discrepancy loss down to similarly low levels. We then give the ē = 0.5 instance another 4500 iterations (purple after red curve) and it could outperform the غ = 0.25 instance. Finally, tried using fully noisy initialization of Ü and the result is much worse than other instances.
In our first benchmark, Taylor vortices Pavloy et al. (2011), we further analyze the sensitivity of our method with respect to the initial guess U. To this end, we first compute U via POD and then corrupt U using a random noise bases U' with each element sampled according to the truncated normal distribution with μ = 0, a= and truncated to range [_1, We then use the following initial guess: Retract(U, ÚT)) where Σ, is scaling diagonal matrix such that each column of ÜΣ has no-m equals to some and ē controls the magnitude of random noise. Here multiplying by DD ensures that our noise is divergence-free In Figure 5, we profile the convergence history with غ = 0.01, 0.05, 0.25, 0.5. Although the noise can drastically change the initial discrepancy loss. all four instances can reduce the loss to similar levels after sufficiently many iterations. Our analysis also implies that the POD baseline provides a good initial guess of U, because a fully noisy initialization of U can lead to worse result.
8

---

Under review as a conference paper at ICLR 2023
coupled(6)
coupledt 5) One-Step (oupled(5) () One Sne-Step
Mode
s/p
Loss-POD
Loss-ak-Iter.
Descrpancc
One-Step
0.05/6
$7.39
24.12
Full-Unrolling 0.05/6
57.39
6.37
One-Step
0.01/9
28.23
15.71
Full- Full-Unroling
0.01/9
28.23
5.38
diteration
Figure 6: The convergence history over 3000 iterations of four instances of learning reduced smoke plume rising trajectory. We use two sets of instances: 0.05, p = and 0.01, p For each set, we compare one-step and full-unrolling mode training.
In the recent work Brandstetter al. (2022), authors proposed two training modes for learning neural PDE solver, one-step training and full-unrolling. One-step training cuts off the gradient after single timestep. while the full-unrooe unlolling mode considers the full gradient of Equation 7 over the entire trajectory. We compare the two modes in Figure 6 in terms of trajectory-wiis discrepancy loss, using our second benchmark scenario, rising smoke plume. Both modes can reduce the loss after 3000 iterations, although there is some initial fluctuation in one-step training, while full-unrolling leads to significantly faster convergence. We use the full-unroling mode for all other examples.
Discregancy
Our Method
POD
125
100
100
160
180
0.3 p.50 6.79 8.n
0.5n
6.
6.5mg
0.7m
(a)
(b)
6.78
(C)
Figure 7: (a): The trajectory-wise discrepancy with respect to 0 for our third benchmark. (bc): The initial (b) and final (c) trajectory- -wise discrepancy with respect to 01, 02 for our forth benchmark.
Our third benchmark involves a spherical smoke plume, with initial diameter 1/3 and speed 1.0 lo-
cated in the center of a [0, 1]² domain, moving in varying directions. We assume the direction of mo- tion is parameterized by the angle 0 [0, 2T] sampled from the initial distribution I 2(( 2 2m)) We use velocity field discretized on 64 x 64 rectangular grid with Dirichlet boundary condition (n = 8320). Our training dataset for the POD baseline contains N = trajectories with evenly sampled 0%, 45" 90"00 .... With
500, öt 0.01, c 0.01, p 36, we run our method for 12200 iterations, taking 72 hours to converge. We then test our method on another 24 evenly sampled = 7.5% 22.5" ,0",.,, , which are not covered by the training dataset (some snapshots can be found in Figure 15 of Appendix H). As plotted in Figure 7a, our method reduces the dis- crepancy by 54.65% on average.
Our fourth benchmark extends the third one by involving two smoke plumes, located (0.5, 0.25) and (0.5, 0.75). The directions motion 0,, 0g e [0, マ are sampled from the initial distribution = [[[[[[]]>]] [?) and we set 0.01, p= = 59. Our training dataset for the POD baseline contains N 25 trajectories with 5 evenly sampled 01.2 0%, 72 144", 216 288*. Other parameters are the same as those of our third benchmark. We run our method for 18000 iterations, taking 72 hours to converge (some snapshots can be found in Figure 16 of Appendix H). Afterwards, we test our method on another 25 evenly sampled 01,2 36 108", 180 3522,3241 that are not covered by the training dataset. As plotted in Figure 7bc, our method reduces the diserepancy by 59. 59.28% on average.
6 CONCLUSION
We propose a model-based approach to fine-tune reduced fluid dynamic systems. Our main idea is to rely the differentiable structure between the state transfer function and the linear subspace bases to minimize the expected trajectory-wiis discrepancy loss, over a distribution of initial con- ditions. By evaluating several simulation benchmarks, we show that method outperforms the POD baseline. As our major limitation, our trajectory prediction has sequential dependence and cannot exploit GPU parallelization. Even with tensor precomputation technique, the training

---

Under review as a conference paper at ICLR 2023
still takes hours a desktop machine, which is orders of magnitude slower than the simple POD or DMD method. In addition, our method uses a linear subspace with limited expressivity as compared with universal neural networks Wu et al. (2021); Hasegawa al. (2020); Lee et al. (2021) used by non-intrusive model reduction techniques. We speculate that using neural networks to represent the reduced hases (/ is possible as done in Li al. (2017), although the orthogonal and divergence-free constraints will be more difficult to enforce, Enforcing these constraints exactly as in Mohan et al. (2020) would compromise the efficacy of reduced time integration.
REFERENCES
Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. Advances in neural information processing systems, 32, 2019.
V.C. Aitken and H.M. Schwartz. On the exponential stability of discrete-time systems with applica- tions in observer design. IEEE Transactions on Automatic Control, 39(9): 1959-192,,,, 1994. doi: 10. 10.1109/9.317135
Ekhi Ajuria Illrramendi, Antonio Alguacil, Michaël Bauerheim, Antony Misdariis, Benedicte Cuenot, and Emmanuel Benazera. Towards an hybrid computational strategy based on deep learn- ing for incompressible flows. AIAA AVIATION 2020 FORUM, pp. 3058, 2020.
Giancarlo Alfonsi. deynolds-averagee stvier equations for turbulence modeling. Applied Mechanics Reviews, 62(4), 2009.
Brandon Amos and Zico Kolter. Optnet: Differentiable optimization as a layer in neural networks. In International Conference on Machine Learning, pp. 136-145. PMLR, 2017.
David Amsallem and Jan Nordstrom. Energy stable model reduction of neurons by nonnegative discrete empirical interpolation. SIAM Journal on Scientific Computing, 38(2):8297-B322 2016.
Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al. Interaction networks for learning about objects, relations and physies. Advances in neural information processing systems, 29, 2016.
Oktay Baysal and Mohamed E Eleshaky. Aerodynamic design optimization using sensitivity analy- sis and computational fluid dynamics. AIAA journal, 30(3)711--72,, 1992.
Gary Becigneul and Octavian-Eugen Ganea. Riemannian adaptive optimization methods. In Inter- national Conference on Learning Representations, 2019.
Thomas Bendokat, Ralf Zimmermann, and P-A Absil. A grassmann manifold handbook: Basic
geometry and computational aspects. arXiv preprint arXiv:2011.13699 2020.
Gal Berkooz, Philip Holmes, and John Lumley. The proper orthogonal decomposition in the analysis of turbulent flows. Annual review of fluid mechanics, 25(1):539 575, 1993.
Lorenz T Biegler, Omar Ghattas, Matthias Heinkenschloss, and Bart van Bloemen Waanders. Large- scale pde-constrained optimization: an introduction. In Large-Scale -DE-Constrained Optimiza- tion, pp. 3-13. Springer, 2003.
Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. arXiv preprint arXiv:2202.03376 2022.
Robert Bridson. Fluid simulation for computer graphics. AK Peters/CRC Press, 2015.
Robert Bridson and Christopher Batty. Computational physics in film. Science, 130(6012):1756- 1757, 2010.
Saifon Chaturantabut and Danny C Sorensen. Nonlinear model reduction via discrete empirical interpolation. SIAM Journal on Scientific Computing, 32(5):2737-2764, 2010.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
10

---

Under review as a conference paper at ICLR 2023
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXiv preprint arXiv:2003.0400 2020.
Thomas Daniel, Fabien Casenave, Nissrine Akkari, and David Ryckelynck. Model order reduction assisted by deep neural networks (rom-net). Advanced Modeling and Simulation in Engineering Sciences, 7(1):1-2) 2020.
Leszek Demkowicz, Tinsely Oden, Waldemar Rachowicz, and Oliver Hardy. Toward universal hp adaptive finite element strategy, part 1. constrained approximation and data structure. Computer Methods in Applied Mechanics and Engineering. 77(1-2):79-11 -122211:91 1989.
Nicola Demo, Marco Tezzele, and Gianluigi Rozza. Pydmd: Python dynamic mode decomposition. Journal of Open Source Software, 3(22):530, 2018.
Tao Du, Kui Wu, Andrew Spielberg, Wojciech Matusik, Bo Zhu, and Eftychios Sifakis. Functional optimization of fluidic devices with differentiable stokes flow. ACM Trans. Graph. 39(6), nov 2020. ISSN 0730-0301.
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of machine learning research, 12(7), 2011.
M. Duponcheel, P. Orlandi, and G. Winckelmans. Time-reversibility of the euler equations as a benchmark for energy conserving schemes. Journal of Computational Physics, 227(19):8736 ISSN 0021-9991.
8752, 2008.
Alan Edelman, Tomás A Arias, and Steven T Smith. The geometry of algorithms with orthogonality constraints. SIAM journal Matrix Analysis and Applications, 20(2):003-333, 1998.
Mare Finzi, Ke Alexander Wang, and Andrew G Wilson. Simplifying hamiltonian and lagrangian neural networks via explicit constraints. Advances in neural information processing systems, 33:
13880-13889, 2020.
Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019.
Gugercin, C Beattie, and AC Antoulas. Rational krylov methods for optimal h2 model reduction.
submitted for publication, 2006.
Kazuto Hasegawa, Kai Fukami, Takaaki Murata, and Koji Fukagata. Cnn-lstm based reduced der modeling of two-dimensional unsteady flows around a circular cylinder at different reynolds numbers. Fluid Dynamics Research, 52(6))06550,, 2020.
Roland Herzog and Karl Kunisch. Algorithms for pde-constrained optimization. Mitteilungen, 33(2):163-176, 2010.
GAMM-
Philipp Holl, Nils Thuerey, and Vladlen Koltun. Learning to control pdes with differentiable
physics. In International Conference Learning Representations, 2020. URL https: / //openzeview. net/
William W Hsieh and Benyang Tang. Applying neural network models to prediction and data anal- ysis in meteorology and oceanography. Bulletin of the American Meteorological Society, 79(9): 1855-1870, 1998.
Byungsoo Kim, Vinicius Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara Solenthaler. Deep fluids: A generative network for parameterized fluid simulations. In Computer graphics forum, volume 38, pp. 50-70. Wiley Online Library, 2019.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv: 1et.. 2014.
SiHun Lee, Kijoo Jang, Haeseong Cho, Haedong Kim, and SangJoon Shin. Parametrie non-intrusive model order reduction for flow-fields using unsupervised machine learning. Computer Methods in Applied Mechanics and Engineering, 384: 13999, 2021.

---

Under review as a conference paper at ICLR 2023
Qianxiao Li, Felix Dietrich, Erik M Bollt, and Ioannis G Kevrekidis. Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(10):103111, 2017.
Yifei Li, Tao Du, Sangeetha Grama Srinivasan, Kui Wu, Bo Zhu, Eftychios Sifakis, and Wojciech Matusik. Fluidic topology optimization with an anisotropic mixture model. ACM Trans. Graph., nov 2022.
Yunzhu Li, Jiajun Wu, Jun- Yan Zhu, Joshua B Tenenbaum, Antonio Torralba, and Russ Tedrake. Propagation networks for model-based control under partial observation. 2019 International Conference on Robotics and Automation (ICRA), pp. 1205-1211. IEEE, 2019.
Julia Ling, Reese Jones, and Jeremy Templeton. Machine learning strategies for systems with in- variance properties. Journal of Computational Physics, 318 318:22-35, 2016.
Beibei Liu, Gemma Mason, Julian Hodgson, Yiying Tong, and Mathieu Desbrun. Model-reduced variational fluid simulation. ACM Trans. Graph., 34(6), oct 2015. ISSN 0730-0301.
Pablo Márquez-Neila, Mathieu Salzmann, and Pascal Fua. Imposing hard constraints on deep net- works: Promises and limitations. arXiv preprint arXiv: 1706.0270000223 2017.
Romit Maulik, Vishwas Rao, Sandeep Madireddy, Bethany Lusch, and Prasanna Balaprakash. Using recurrent neural networks for nonlinear component computation in advection-dominated reduced- order models. arXiv preprint arXiv:1909.0914 0444 2019.
Aruna Mohan and Glenn H Fredrickson. Reduced-order computational model for the molecular dynamics simulation of entangled polymers. arXiv preprint arXiv:2009.00211 2020.
Arvind T Mohan, Nicholas Lubbers, Daniel Liveseu, and Michael Chertkov. Embedding hard physical constraints in neural network coarse-graine goaseegraining of 3d turbulence.
arXiv preprint
arXiv:2002.00021, 2020.
Parviz Moin and Krishnan Mahesh. Direct numerical simulation: a tool in turbulence research. Annual review of fluid mechanics, 30(1):539-578 1998.
Alexander Morozov, Davide Zgyatti, and Petr Popov. Equidistant and uniform data augmentation for 3d objects. IEEE Access, 10:3766-3774, 2021.
Richard M Murray, Zexiang Li, and S Shankar Sastry. A mathematical introduction to robotic manipulation. CRC press, 2017.
Sina Ober-Blöbaum and Christian Offen. Variational learning of euler-lagrange dynamics from data. arXiv preprint 12Xiv:2112.12619, 2022.
Dmitry Pavlov, Patrick Mullen, Yiying Tong, Eva Kanso, Jerrold E Marsden, and Mathieu Desbrun. senccture-preserving discretization of incompressible fluids. Physica D: Nonlinear Phenomena, 24006))443-458, 2011.
Barak A Pearlmutter. Learning state space trajectories recurrent neural networks. Neural Com- putation, 1(2):263-669, 1989.
Titus Petrila and Damian Trif. Basics of fluid mechanics and introduction to computational fluid dynamicy, volume 3. Springer Science & Business Media, 2004.
Lukas Prantl, Boris Bonev, and Nils Thuerey. Generating liquid simulations with deformation- aware neural networks. In International Conference on Learning Representations, 2019. URL
.et/ forum widddid=-Hy@ss jON.GBBjJ30 9F
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informe neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
12

---

Under review as a conference paper at ICLR 2023
Börner Ralph-Uwe, Oliver G Ernst, and Klaus Spitzer. Fast 3-d simulation of transient electro- magnetic fields by model reduction in the frequency domain using krylov subspace projection. Geophysical Journal International, 173(3):766-780, 2008.
Clarence W Rowley. Model reduction for fluids, using balanced proper orthogonal decomposition. International Journal of Bifurreation and Chaos, 15(099)997-1113, 2005.
Clarence W Rowley and Scott TM Dawson. Model reduction for flow analysis and control. Annu. 2017.
Rev. Fluid Mech,
Clarence W Rowley and Jerrold Marsden. Variational integrators for degenerate lagrangians, with application to point vortices. In Proceedings of the 41st IEEE Conference on Decision and Control, 2002., volume pp. 1521-1527. IEEE, 2002.
Rubens Sampaio and Christian Soize. Remarks on the efficiency of pod for model reduction in non-linear dynamics of continuous clastic systems. International Journal for numerical methods in Engineering, 72(1) ((22):22443 45, 2007.
Connor Schenck and Dieter Fox. Spnets: Differentiable fluid dynamics for deep neural networks. In Conference on Robot Learning. pp. 317-335, PMLR, 2018.
Peter J Schmid. Dynamic mode decomposition of numerical and experimental data. Journal of fluid mechanics, 656:5-28, 2010.
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339-1364 2018.
Georg Still. Lectures on parametric optimization: An introduction. Optimization Online, 2018.
Kunihiko Taira, Steven L Brunton, Scott TM Dawson, Clarence W Rowley, Tim Colonius, Bev- erley J McKeon, Oliver T Schmidt, Stanislav Gordeyev, Vassilios Theofilis, and Lawrence Ukeiley. Modal analysis of fluid flows: An overview. Aiaa Journal, 55(12):401344441 2017.
Tetsuya Takahashi, Junbang Liang, Yi-Ling Qiao, and Ming C Lin. Differentiable fluids with solid coupling for learning and control. In Proceedings of the AAAI Conference on Artificial Intelli- gence, volume 35, pp. 6138-6146, 2021.
Qingyang Tan, Zherong Pan, Breannan Smith, Takaaki Shiratori, and Dinesh Manocha. N-penetrate: Active learning of neural collision handler for complex 3d mesh deformations. In International Conference on Machine Learning, pp. 21037-21049, PMLR, 2022.
Kiwon Um, Robert Brand, Yun Raymond Fei, Philipp Holl, and Nils Thuerey, Solver-in-the-lopp: Learning from differentiable physics to interact with iterative pde-solvers. Advances in Neural Information Processing Systems, 33:611 111111661210 6122, 2020.
Zheng Wang, Dunhui Xiao, Fangxin Fang, Rajesh Govindan, Christopher C Pain, and Yike Guo. Model identification of reduced order fluid dynamics systems using deep learning, International Journal for Numerical Methods Fluids, 86(4):255 2268 2018.
Pin Wu, Siquan Gong, Kaikai Pan, Feng Qiu, Weibing Feng, and Christopher Pain. Reduced order model using convolutional auto-encoder with self-attention. Physics of Fluids, 33(7):077107, 2021.
Ping Yang, Yao-Lin Jiang, and Kang-Li Xu. A trust-region method for h2 model reduction of bilinear systems on the stiefel manifold. Journal of the Franklin Institute, 356(4):2258 -273,, 2019. ISSN 0016-0032.
Thomas Zang and Lawrence Green. Multidisciplinary design optimization technique--imllicaionns and opportunities for fluid dynamics research. In 30th Fluid Dynamics Conference, pp. 3798, 1999.
Jun Zhang. Fast and high accuracy multigrid solution of the three dimensional poisson equation. Journal of Computational Physics, 143(2):449-44) 4611 1998.
Weiwei Zhang, Jiaqing Kou, and Ziyi Wang. Nonlinear aerodynamic reduced-order model for limit- cycle oscillation and flutter. Aiaa Journal, 54(10)33304-3311 2016.
13

---

Under review as a conference paper at ICLR 2023
A DISCRETE ENERGY PRESERVATION
We prove that energy preservation and time reversibility hold in a time-discrete setting. Lemma A.1. The tensor Cksj antisymmetric.
Proof. This follows from the definition of Ckaz:
Cuis = SM JU (UK, VxU,XX
-VU:- U (VU-- VU
=- --JM UJ -UTUU-V SM (U,, V V = -Gjkk
where we used elementary vector identity that (V A)x B = B.
Using the antisymmetry of Cuig, we can show that trapezoidal rule is indeed energy preserving. Lemma A.2. For any 2" satisfying the trapezoidal rule, ||=*| |=|
Proof. Multiplying the lefthand side of Equation 4 by # zy. and summing over k. we have:
x+ 11
|z||2
Σ
Ckij
2k
2
2k
2.
+
24
23
+
öt
2
2
2
kay
122 |0||a
+
Σ
(Ckis + Cjass)
+
2.
+ 2. 23
+
2j
|= -1 || 0,
ôt
2
2
2
ôt
kdy
from which result follows.
Next, we show that the trapezoidal integrator (Equation 4) must have solution by a proper choice of sufficiently small 8t.
Lemma A.3. There exists a sufficiently small 8t such that Equation be solved for z" via the following negative gradient flow:
with initial guess z²
Proof. We consider the Lyapunov candidate V(=*) +1{==)] on the ball B,(x) f=+11=-- -
r). The negative gradient flow satisfies:
v(z") =
s- (1)--))((==+((--=))))((++++)++++
Now since the eigenvalue of a Hermitian matrix is a Lipschitz function of matrix entries Golub & Van Loan (2013), we must have:
e(I=IIr) & <(VC(==)vC(==))) s ((=1,rr)
for some p, and any =t e (z). Combining the above estimation, we have:
(1 (0-000000000
Obviously, with sufficiently small öt, we have x for some (0, and B, (2). Next, consider the boundary case 2" E EB,,(2), where we have:
V(**) - + |
StC
25t [IC(=*))
(11-)t)²* (ot² 00) ((=)) 13
and we can choose sufficiently small öt such that V(=*) V(2) for all zt e OB, (z). Our result follows from the exponential stability condition Murray et al. (2017).
14

---

Under review as a conference paper at ICLR 2023
In practice, however, continuous gradient flow cannot be realized, but similar argument as Lemma A.3 can be used to show that the Newton Newton-Raphson method is guaranteed to converge when minimizing V(=*) under sufficiently small öt:
Lemma A.4. There exists sufficiently small öt such that Equation 4 be solved for 2" via the Newton-Raphsse method:
( =2(-1) (aneld-))--((((aa1
with initial guess z 20) = 2. Here we use superscript with bracket to denote iteration index.
Proof. Consider the reduction of Lyapunov candidate V(2)
after one iteration, have:
ôt 1(2(1-1)T (2
V(z(1) =I - f(z (1=
Σ k
Hx Hx( (Vf(2( Ck. O (((2(111)))1
By similar argument as in Lemma A.3, we can choose sufficiently small öt such that:
((Hx(2((11)) s +() V(2(0) pôt² p ((1=11,)) 1f(2((-1))
as long as 2 (d-1) (z)... We can also choose sufficiently small 5t such that:
V(z(0) ()(((())) (2(d-1) I[) (d-1) =1,
(8)
for some e (0, 1). Next, we consider the Hessian of Vz²))
v²V(=")
8
where R(²*) is smooth. symmetric matrix function. We can further choose sufficiently small öt such that (x=*) is 1/2-strongly convex and for any 2t Byr (z):
V(=* ()=) (+ -z)
By the smallness of öt, we have:
v(=t)> V()) (V
Vzt E B3r (z) (2(d-1) E (x)=) A (2 min(1,
(9)
Combining Equation 8 and Equation 9, we have for small enough öt:
(2 Bx(z) (V()(4) seV(2(a-1) (2(1-) (()()))((((--)) mmin(1,2/2))
Our result follows by choosing sufficiently small öt such that V(x(0)) S min(1,,2/22 and invoke the discrete exponential stability condition Aitken & Schwartz (1994).
Note the choice of öt is only dependent |2|| and r, which can be used to show that the timestep size can be fixed throughout the trajectory for time reversible fluid systems:
Corollary A.S. Given an initial condition zo, an energy preserving discrete trajectory can be com- puted by repeatedly solving Equation for zx using a fixed timestep size ől via the Newton Raphhoo method.
Proof. This result can be derived by induction on two facts: | = by Lemma A.2; 2) To solve for 2 St can be determined function 50(1|===11,,) II,)) by Lemma A.4.
n
15

---

Under review as a conference paper at ICLR 2023
B DISCRETE TIME REVERSIBILITY
The above result guarantees energy preservation throughout the trajectory. We now move on to show time reversibility the discrete setting:
Lemma B.1. There exists a sufficiently small öt, such that for any = & B, (0), the negative gradient flow Equation 4 defines invertible map from to z².
Proof. Following the same argument as in Lemma A.4, we can choose sufficiently small öt such that V(=") is strongly convex when restricted to (2) and the map = (=) argminV(2") is
well-defined and differentiable Still (2018). The derivative of function z (z) can then be derived via the implicit function theorem as:
V2" (2) (1 [1-
By Lemma A.2, we know that x- e B, (0) as well. Again by the lipschitz continuity of singular values, we can choose sufficiently small öt such that dct(Vz" (2)) throughout (,) and our result follows by the inverse function theorem.
Lemma B.1 can also be extended to the entire trajectory via induction:
Corollary B.2. Given an initial condition 2 e B,(0) for some r, an energy preserving discrete trajectory can computed by repeated solving Equation 4 for zh using a fixed timestep size öt, such that the resulting map 2 (=") is invertible.
Proof. By induction on Lemma A.2 and Lemma B. we know that 2* (*(a11) is invertible for any k > 0 and our result follows by composition of invertible functions.
Ou Meliol
JD
Figure 8: illustrate four man- ifolds: St(n,p) for the velocity bases; Si(m.,) for the divergence-free veloe- ity bases; Gr(n.p) for the velocity sub- space; Gr(m,
the divergence-fres velocity subspace. Our method main- tains e St( St(n,)) p) and represents the gra- dient as some VoL St (SSnn,)),,,,)))))))))p() p). which is both memory efficient and computa- tionally tractable.
pi Ca(mmap)
C DERIVATIVE FORMULATION
In this section, analyze the differentiability of lifted transfer function Equation To com- pute derivatives of the forward dynamic function with respect to the bases U, we need to utilize the implicit function theorem and special representation of the bases as manifold point, which cannot be exploited by automatic differentiation. First, we show that the function is well-defined on the manifold Gr(m,p) via the following lemma:
Lemma C.1. The listed transfer function Equation can be written as a fanction of (ū, pm).
Proof. By the incompressibility of bases Ú = T we have: P DDDm) D". Plugging this into Equation and we have the follow rewrite:
87 (ū, pm)): {
Piu+= P
CCP pett))
from which our result follows. We can derive the original definition (Equation 5) by multiplying the second equation by T from the left.
n
16

---

Under review as a conference paper at ICLR 2023
Although the function is well-defined, the complexity of its derivative computation relies an efficient representation of bases. A straightforward representation is to use matrix pm and consider the function 8* (U,Pm) However, this representation requires storing the large matrix pm which is computationally impractical. In this section, we exploit equivalent manifold representations to derive the computationally tractable formulas for the derivatives of arbitrary loss functions L a i²- The relevant manifolds are illustrated in Figure 8, We first derive the partial derivative 0o" /0r via the implicit function theorem:
Oo Oe
[ULL---TCIe"" P..
(10)
The inverse of the system matrix above is well-defined when the timestep size is sufficiently small according to Appendix A. It can be verified that the above derivative is invariant to the orthogonal basis transform. Next, we derive the partial derivative with respect to pm E GG(mm,)) p). We denote UTTh as the complement of Umh and QT = (Um, UT E O(m). An element of Tpm Gr(m, p) can be identified with matrix dB E ((m--)<x via:
dBT dB
dpm = Q"
We can lift Gr(m, to St(m,,) via the map Ts((mmp)--Cr(mmpp (1 = Umm [ūm]". Under this map, an element dt'm E St(m.)) p) horizontal of TT.... Gr (p) must satisfy the condition dūmm UTB (we refer readers to Bendokat et al. (2020) for the derivation). Representing gradient as some dB's is the most memory efficient method, since the dimension of Gr(m, p) equals that of dB. However, we have to multiply dB with U!" and then with D to recover divergence-free velocity bases, while computing either U or D is intractable. Instead, we choose to work with dŪ directly and rely on the following result that establishes a connection between dB and dU:
Lemma C.2. For divergence-free velocity bases U, a direction dŪ belongs the tangent plane of ((, p) n S(n,p) at U if and only if dÙ E D(n,p) and U dU O.
Proof. If dū belongs to the tangent plane, then it must satisfy dū = "dD for some dtnh = UTM dB, so dÙ' E p).,p). Further, t' = 0 DŪT dB = [ 1 UTH dB = O. Conversely, di D(n. p) implies dū = Dd for some dim Further, ÜT du = 0 implies [tm] T dū = 0, which in turn implies dL = CTT dB for some dB.
Suppose we have loss function L Tz (ū. Pmm with i as the constant, we can composite the loss function with the map function is the intersection of D(n.p) p) and p). which is an embedded sub-manifold of RT*P. In order to calculate the gradient the manifold, we can smoothly extend the composite function to the entire Ruxp, calculate the Fuctidean-spac. gradient denoted by G e Rup,, and then project the gradient onto the tangent space. Such projection is defined by Lemma C.2 as:
(U) DTOUTD pm. The domain of this composite
VeL=P,DD"G.
(11)
where multiplying by DDT ensures VoL E D(n,p) and multiplying by P. ensures VTV = O. Note that, although computing the entire 1) is intractable, evaluating DDTG G is tractable. Indeed, this involves projecting each column of G into the divergence-free vector subspace, which can be calculated by solving discrete Poisson's equation Petrila & Trif (2004) via a sparse linear solve at a complexity of Onn)) Zhang (1998), where w 2 1 depends on the numerical linear system solver. Therefore, the entire projection has cost of Onn'p), puting 1 being O m)) We refer readers to Appendix C.I for the derivation of Euclidean space gradient G. The computation of v. over long trajectory with T >> p timesteps is rather efficient. Indeed, we can precompute and accumulate G for each timestep, and finally apply divergence-fire -ree projection operator to compute V C.C, the total cost of which is O(n p+ + Tnp + Tp²).
compared with the complexity of com-
C.1 DERIVATIVE FORMULATION IN EUCLIDEAN SPACE
We derive the formula for G in the following lemma:
17

---

Under review as a conference paper at ICLR 2023
Lemma C.3. If introduce the third order tensor:
Pa Σ U2, 2/ + 23 +
EC( 2 2,
2
2j
[C(Ua.Ur.e)) U,. 2 + 2
and consider an arbitrary differentiable function [(ū), then the Euclidean space gradient G of function L o o- (v,U) with respect to U defined as:
= - ( v ==]" - -V
(12)
Proof. Assuming 0 is fixed, we first derive some useful fundamental results:
d= = =[[ D' ú = dB [ŪT]] D'
d PT) =d DU"> DdU 2 DU
D UT dB [mm]" +Ü"d [UT T D
dBT
= =πQm
[ DJ ū= dPr = -dP
dB
Plugging & into the first-order expansion Equation and we have:
[I StVC(= 4 dū StvC(=*)] [T- StTC(=
where: denotes tensor contraction of the last two indices. The remaining derivation follows the chain rule:
dū =Üd=* + d0 (tTT=++ (ura=
+ dú]+
dÙ |= -
dC =v do tr(dU G).
By comparing the two sides of the last equation, result follows.
C.2 ALTERNATIVE LIFTED FUNCTION
The above derivation is based on the definition of 5* (ū, U) in Equation 5, which assumes that the orthogonal component is kept across timesteps. An useful alternative is assume that the orthogonal component is discarded, which is:
5* o'(,,0) 0==))===(000)))
(13)
By a similar argument, we can derive the following derivatives for Equation 13:
Ov
De [
V
D DECOUPLED REDUCED-ORDER MODEL.
We observe that energy preservation and time-reversibility discussed in Appendix A only requires the tensor Ckiy to be antisymmetric. other words, the construction of the tensor Cat via Equa- tion 3 is not necessary. We speculate that using learned antisymmetric tensor Ckij can expose larger search space, leading to better match with the full-order model. We denote such model as decoupled reduced-order model, where Ckij are separate decision variables not constructed from
18

---

Under review as a conference paper at ICLR 2023
Benchmark
coupled #variable (np)
decoupled #variable (mp p²)
Taylor Vortices
8192
8/1 1/11/1
05536/90122/131022000
66048/914433133316622200 48600/73305/ 24335522240
Plume Rise
8064
6/9/15/26
483844/25566/ 20960/06644
Plume Obstallss
7416
5/8/16/30
37080/59328/1 18666/222480
Spherical Plume
8064
290304
Two Plume
8064
475776
summarized number of decision variables each example. the coupled case, our decision parameter is U having up variables. In the decoupled case, our decision variables are U, Ckis having np+ + p² variables.
Table 2:
U. The formula for 0o /00 Equation 10 stays the same and the formula for G takes the following simpler form:
G -==0"""01
=1" -oL
Finally, the derivative with respect to Cisk reads:
= [I OVVC )
aC
L 2|t,, 5it;
Laj
2
where we have projected the derivative onto the antisymmetric subspace. On the downside, there universally valid öt to make our objective function globally differentiable for all Ü and Cıajy because discrete time reversibility requires sufficiently small öt that depends on Caiy- Empirically, however, have not observed any convergence issue. In Figure 9, we compare the coupled and de- coupled versions on the Taylor vortices and the smoke plume benchmark, their convergence histories almost identical. Therefore, recommend always using the coupled model due its theoret- ical differentiability guarantee. In Table 2, we summarize the number of decision variables in our various experiments.
coupled(11) decoupled(11)
coupled(9) decoupled(s)
Discrpancyy
Benchmark
Loss-Init.
Loss: 10k- Iter.
Taylor Vortices Coupled
3.93
0.39
25
(b)
Taylor Vortices Decoupled
3.93
0.42
Plume Rise Coupled
28.22 28.23
5.38
Plume Rise Decoupled
28.23
3.84
Elapsed Time eni
Elapsed Time thr)
Figure 9: We compare the performance of coupled and decoupled versions on the Taylor vortices benchmark (a), with = 0.01 and p= 11, and the smoke plume benchmark (b), with = 0.01 and p
E COMPARISON WITH ALTERNATIVE Loss
To highlight the effectiveness of our physics correctness loss, we conduct a comparison with two other loss functions: the Li and Le losses defined as:
L1 Li 0) 10* -6/1 (2(0*, 0) 1 u||2,
where we denote è as the velocity generated by the groundtruth fullspace fluid simula- tor Pavlov et al. (2011). We note that these loss functions are impractical for large-scale test cases because they require solving for the groundtruth data of a different initial condition dur- ing each iteration of training. Therefore, we choose to only evaluate them on our first three benchmarks in Table 1, where there is only single trajectory so i can be precomputed. For these benchmarks, we both train and evaluate them on the three losses Layn: L1, L2, and sum- marize the results in Table 3. We also plot the convergence history of the first benchmark (Taylor Vertices) in Figure 10. Our plots show that, when the first benchmark is trained using L1,2, L1,2 will both decrease by at most 64%, but our Layn can increase drastically by at most 1083%. Instead, when trained using Layas L1 will increase or decrease by at most 3.3% but our Laya can decrease significantly by 76.8%. Considering these properties and the fact that L,, Lg is impractical to compute by requiring the groundtruth data, we conclude that our Layn is overall more practical in training reduced fluid systems.
19

---

Under review as a conference paper at ICLR 2023
Taylor Vertices (p=8)
Plune Rise (p=b)
Plune Rise+Obstacle (p=5)
Train
Bval.
Cays
(10 Co(10 ")
Train
Eval.
Liyn
L1 C, +)
Train
Eval.
Lap
( 4) C2(10 +)
Las
0.58
1,23
0.27
6.37
2.37
0.79
Cam
10.48
1.92
0.53
C.
37.81
0.44
0.04
C
$7.50
0.67
0.07
Ci
132.68
0.49
0.06
L
165.41
0.27
0.06
42
114.20
0.43
0.05
Ca
87.83
0.69
0.06
Table 3: We evaluate our first three benchmarks when trained and evaluated using Layn, Li Lz.
Ic
1e-3
2.5
Our Loss
1.2
Our Loss
1. Loss
Our Loss
Loss Loss
6 Loss
5,2 Loss
0.8
200
Our
1a Loss
0.6-
1.0
50
10
0.5
0.00
0.0
0.5
1.0
1.5
(a)
0.5
1.0
1.5 led
(b)
0.3
1.0
1.5
2.0
(c)
Aiteration
1n4
#larraiinn
Alienation
Figure 10: For our first benchmark (Taylor Vortices), we plot the convergence history when trained using Layn (a), Li (b), and L2 (c). The scale of Cayn shown on the left and Lı, shown on right of each plot.
F COMPARISON WITH DMD
We have shown that our method works best with POD initialization. In this section, we con- duet additional experiments with DMD. DMD extends POD by assuming that the data is gen- erated from a linear dynamic system. DMD can be used both as an intrusive and non-intrusive method. In the intrusive mode, we use DMD to compute a bases Ü and compute Cki from U via Equation namic system as the surrogate. To evaluate the performance of DMD, we use two metrics. For the intrusive DMD,
In the non-intrusive mode, we simply use the DMD-assumed' linear dy-
use our physics correctness loss Equation 6. Unfortunately, our physies correctness loss is not suitable for evaluating non-intrusin methods that can be non-reversible. Indeed, is always possible to let Layn by setting ō" O. Therefore, we also measure as an indication dynamic system stability.
the energy gain Ae (1771 1/ 77
We perform the experiments using the open source DMD library Demo et al. (2018) on our first three benchmarks. Their results are shown Table 4. The results show that the per- formance of intrusive DMD is worse than either POD or our method, in terms of the physics correctness loss. This because the main assumption of DMD, i.e., the dynamic system being linear, is invalid for the bilinear dynamic system Equation 3. Instead, POD does not make any assumption on the time dependency between frames and serves as a better initialization for our method. On the other hand, the non-intrusive DMD leads to better performance in terms of Layn but the dynamic system tends to be rather unstable due to a drastic energy gain of 1.9 -73.xx.
Loya/ Ae
Benchmark
P
POD
I-DMD
NI-DMD
Ours
Taylor Vortices
8
4.84/0
20.81/0
4.77/4.77
0.58/0
Plume Rise
6
57.04/0
127.18/0
2.21/1.94
6,37/0
Plume Rise+Obstacle
133.16/0
245.20/0
2.32/73.3
10.48/0
Table 4: We compare the POD baseline and our method with intrusive DMD (I-DMD) and non-intrusive DMD (NI-DMD) in terms of trajectory-wise physics correctness loss and energy gain.
G COMPARISON WITH PINNS
We conduct comparisons with PINNs Raissi et al. (2019). PINNs was originally designed for solving PDEs, while our divergence-free Navier-Stokes equation is an DAE. In order to ex- tend PINNs to handle DAE, we learn neural network DAE solution function, denoted as
20

---

Under review as a conference paper at ICLR 2023
g,(x,g,t) (Ua, Bys A) and represented MLP with 3 hidden layers each having H neurons and Tanh activation function, and minimize the following physics violation loss:
| V x VAI vl
We also enforce additional temporal and spatial boundary conditions loss functions. All the loss functions have weights equal to 1. For fairness of comparison, we use the same training data for both our method and PINNs. Note that our method uses grid-based spatial discretiza- tion, so we use all the grid centers as spatial samples of training data and we sample the tem- poral domain predict trajectory of the same length as our method, i.e. Tôt. We use Adam as our optimizer and we train both methods on CPU for 24 hours. Since PINNS can lead non-divergent-free velocity fields, we measure the accuracy of both methods via three metries: Layns Ae, and av- erage divergence error: 0 - 0" oooo where ü* is the closest divergence-fre velocity field c. The results summarized in Table 5.
regular interval of öt 0.01, which equals our timestep size. We aim to
PINNs(H = 64)
PINNs(H = 128)
Ours
Benchmark
Layn/De De
| - v" 80
Layn/ Ae 8.96/0.64 64
|| u 11
Layn/ De 0.58/0
Taylor Vortices
17.46/0.32
0.000027 0.005159
0.000018 0.004466
Plume Rise
6.74/1.06
6.29/0.91
6.37/0
0
Table 5: We compare our method and PINNs in terms of Layn and 0 0² I.0
PINNs mostly perform worse than our method in terms of Ldyn- In the Taylor Vortices bench- mark using H 128, the Laya metric generated by PINNs is slightly better than our method. But this again because Layn is only designed for measuring time-reversible flows, which is not effective metric for comparing reversible and non-reversible flows due to its trivial so- lutions. Such trivial solutions are indeed exhibited in PINNs, as illustrated in Figure 11. After very short period of time, the solution predicted by PINNs become significantly smeared out and meaningless.
(A)
(b)
(6)
Figure 1I: We compare frames generated by groundtruth (a), PINNs(HI = 128) (b) and our method 0.0001, 25) (c) the Taylor Vortices benchmark. After very short time period, the results generated by PINNs become significantly smeared out and meaningless.
H ADDITIONAL RESULTS
We demonstrate additional experimental results. Some snapshots of our 4 benchmark scenarios shown in Figure 12, 13, 14, 15, and 16, respectively.
21

---

Under review as a conference paper at ICLR 2023
(b)
Figure 12: Velocity magnitude field snapshots of the Taylor vortices benchmark, generated by full-order model (a), our method with 0.0001 and 25 (b), and POD with 0.0001 and 25 (c).
coupled(6) coupled(9)
Method
Discrppancy
screpancc
coupled(15)
couplad(26)
#Bases
(b)
Elapsed Time (hr) (c)
Figure 13: (a) Velocity magnitude field snapshots of the smoke plume benchmark, generated by full-order model (top row), our method with 0.0001 and p 26 (middle row), and POD with 0.0001 and 26 (bottom). (b) Trajectory wise discrepancy loss of POD and method, under different p. (c) The convergence history of method over 24 hours.
22

---

Under review as a conference paper at ICLR 2023
Figure 14: Velocity magnitude field snapshots of the smoke plume benchmark with spherical obstacle,
generated by full-order model (a), our method with 0.0001 and p 26 (b), and POD with 0.0001 and p= 26 (c).
Figure 15: Velocity magnitude field snapshots of the spherical plume benchmark, generated by full-order model (a), our method with 0.01 p= 36 (b), and POD with 0.01 p = (c). The plume moves along 7.5" (arrow), which is not covered by our training dataset.
23

---

Under review as a conference paper at ICLR 2023
Figure 16: Fluid dye field snapshots of two cubical plume henchmark generated by our method (time flows to the right). The fluid dye field is initialized as the two cubes (black) and passively advected by the velocity field. With 0.01 and only 59 bases, we can predict family of trajectories with the two plumes moving different directions (arrow).
24


Reviews:
[{"review_id": "sxMcpymXGv", "reviewer": "Reviewer_eGGS", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a method to learn the time-reversible, nonlinear dynamics of an incompressible fluid. For computational tractability, the authors learn a reduced model over a limited set of initial conditions with tensor precomputation. The authors show that their approach predicts solutions significantly closer to the full order method than conventional methods.", "summary_of_the_review": "The paper is clearly written and well organized, and appears to extend an existing method in a nontrivial way.", "strength_and_weaknesses": "The authors address the challenging problem of predicting the state evolution of high-dimensional, time-reversible, nonlinear fluid dynamic systems; however, I am uncertain how to understand the magnitude of the contribution. While they demonstrate strong performance against POD; as I understand, they also initialize with POD. I am also unsure why the authors choose to evaluate against a nondescript POD baseline rather than DEIM or another intrusive reduced order method noted in the related work. Moreover, as the authors mention, the method takes considerably longer to run than the baseline approach.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "jdrvzijBlc", "reviewer": "Reviewer_J4j9", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper introduces a machine learning approach to identify model-reduced dynamical systems. The main idea is to use stochastic (Riemann) optimization which minimizes the expected trajectory-wise reduction error over a distribution of initial conditions. This is achieved by formulating the fluid dynamics problem as an invertible state transfer function. \nThe paper identifies fluid velocities as a point on a manifold where the reduced trajectory depends on the choice of the subspace of fluid velocities. This map from subspaces to reduced trajectories is shown to be globally differentiable.\nThe paper compares against Proper Orthogonal Decomposition as an example of data-driven sub-manifold model.\n\n", "summary_of_the_review": "The idea and the topic in general are quite relevant.  For the current version of the paper, it is however hard to place this work both wrt to existing numerical methods, but also with respect to neural solvers, or neural surrogates. That makes it really tough to properly weigh the pros and cons wrt to e.g. Proper Orthogonal Decomposition, Dynamic Model Decomposition, but also operator learning methods or PINNs.", "strength_and_weaknesses": "Strengths:\n-\tLearning reduced fluid dynamics is a very relevant research direction\n-\tThe idea is fundamentally interesting, i.e. to use automatic differentiation to optimize numeric algorithms which have pre-calculated parameters that are used as initialization.\n-\tI appreciate the time reversible fluid dynamics section, and the transition to reduced model optimization. \n-\tThe results over the baseline look quite convincing, where it is however hard to place the baseline.\n\nWeaknesses:\nMy general comment is that it is hard to place this work both wrt to existing numerical methods, but also with respect to neural solvers, neural surrogates. That makes it really tough to properly weigh the pros and cons wrt to e.g. Proper Orthogonal Decomposition, Dynamic Model Decomposition, but also operator learning methods or PINNs. Furthermore, I would like to see a test against a model which does not preserve time reversibility. \n-\tThe only comparison is against Proper Orthogonal Decomposition (POD) which as stated in the introduction “is flawed in that it ignores the temporal dependence of state variables”. Naturally the question is for example how Dynamic Model Decomposition (DMD) is doing for the problems.\n-\tAlso from a Deep Learning perspective, e.g. Figure 4 could be learned by an operator learning that takes the initial state as input and outputs the state of the system after e.g. every 100th step. On the other hand, if done with e.g. PINNs this could be performed completely without training data. Just using the boundary conditions for the loss and the equation for the residual loss. It would be really interesting to have a comparison wrt to speed and accuracy. \n-\tPyTorch and more specifically the automatic differentiation capability of Pytorch is used as an optimization tool. It is however hard to judge from the paper what exactly the parameters are that are optimized. How many parameters are optimized for the different problems? Algorithms 1 and 2 in the appendix help a great deal but it took me a great deal to scroll back and forth. Is it possible to place the algorithms in the main paper and write refer to the important parts in Section 3,4?\n-\tThe loss should be a bit more central to the writing of the paper. In Fig 1 the most important components are the reduced-order fluid model and the trajectory-wise discrepancy loss. However, Sections 4.3 takes a huge part of the main paper and is actually hard to follow. In my opinion the readability of Section 4 can be improved by focusing more on the relevant parts and not let the reader figure out what they are. \n-\tWhy is e.g. JAX not used for optimization? Did the authors consider that? \n-\tHow do optimization time vs inference time relate between the different models?\n-\tHow do the individual components of the algorithm relate to the performance, e.g. how does the performance change for other loss terms? \n-\tIs it correct that for the first benchmark a single trajectory is used, whereas in the second benchmark more trajectories are used for optimization? On which trajectories is the testing done afterwards? Shouldn’t there be more than just one trajectories to get a better performance estimate?\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "CLoth-9vrK", "reviewer": "Reviewer_Zt39", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes to refine orthogonal bases, as commonly used for \"proper orthogonal decomposition\" (POD) algorithms with a stochastic, data driven optimization. The paper focuses on a class of reversible, circulation-based flow descriptions. The improvements in terms of reconstruction errors of the optimized basis are evaluated for Taylor vortices and different rising plume variations.\n", "summary_of_the_review": "That being said, I think the direction of the paper is interesting. Improved orthogonal bases could definitely find applications, despite potentially being not up-to-par with learned representations. In its current state, I think the work is not really suitable for ICLR, but potentially other more CFD focused venues. For ICLR, I would expect a stronger focus on learning, and actual demonstrations that the proposed method is up to par with the current state of the art for learned fluid simulations.\n", "strength_and_weaknesses": "As strengths of the paper I see the clear reductions in terms of reconstruction errors compared to regular POD approaches. This is clearly demonstrated, and it's neat to see how much difference this fine tuning step can make. Also, I'm not aware of existing works using this  optimization procedure.\n\nAs weaknesses, I primarily see the narrow scope of the work, and its somewhat detached nature from learning approaches. The work relies on a Eulerian DEC discretization for reversible flows, instead of targeting general flow-like, divergence free motions. These are not too widely used, as far as I can tell. In addition, the focus on improving POD methods further limits the scope of applications. POD methods are popular as baselines, and for special applications, but not really state of the art. In addition, the clear lack of comparisons with learned baselines highlights that the focus on basis constructions is not a big topic at ICLR. From an ICLR paper I would have expected links and evaluations of how this method fares in the context of other learning-based approaches. This is mentioned in the outlook, but unfortunately missing in the current version.\n\nI would also recommend to better justify the focus on reversible flows. As mentioned above they're not widely used, and not very relevant for game applications. After reading about games as targeted areas of applications in the introduction I was expecting the paper to take a different direction.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}]

Author rebuttal:

Decision: Reject