Physics Informed Machine Learning of SPH: Machine Learning Lagrangian Turbulence
Authors: Michael J Woodward
Yifeng Tian
Criston Hyett
Chris Fryer
Daniel Livescu
Misha Stepanov
Michael Chertkov
Venue: ICLR 2022 Submitted (2022)
Keywords: Physics Informed Machine Learning
Smoothed Particle Hydrodynamics
Sensitivity Analysis
Differentiable Programming
Mixed Mode Automatic Differentiation
Deep Learning
Turbulence
Lagrangian Fluid Simulation.

Abstract: Smoothed particle hydrodynamics (SPH) is a mesh-free Lagrangian method for obtaining approximate numerical solutions of the equations of fluid dynamics, which has been widely applied to weakly- and strongly compressible turbulence in astrophysics and engineering applications. We present a learn-able hierarchy of parameterized and "physics-explainable" SPH informed fluid simulators using both physics based parameters and Neural Networks as universal function approximators. Our learning algorithm develops a mixed mode approach, mixing forward and reverse mode automatic differentiation with forward and adjoint based sensitivity analyses to efficiently perform gradient based optimization. We show that our physics informed learning method is capable of: (a) solving inverse problems over the physically interpretable parameter space, as well as over the space of Neural Network parameters; (b) learning Lagrangian statistics of turbulence; (c) combining Lagrangian trajectory based, probabilistic, and Eulerian field based loss functions; and (d) extrapolating beyond training sets into more complex regimes of interest. Furthermore, our hierarchy of models gradually introduces more physical structure, which we show improves interpretability, generalizability (over larger ranges of time scales and Reynolds numbers), preservation of physical symmetries, and requires less training data.

Under review as a conference paper at ICLR 2022
PHYSICS INFORMED MACHINE LEARNING OF SPH: MACHINE LEARNING LAGRANGIAN TURBULENCE
Anonymous authors
Paper under double-blind review
ABSTRACT
Smoothed particle hydrodynamics (SPH) is a mesh-free Lagrangian method for obtaining approximate numerical solutions of the equations of fluid dynamics, which has been widely applied to weakly- and strongly compressible turbulence in astrophysics and engineering applications. We present a learn-able hierarchy of parameterized and "physics-epplainable" SPH informed fluid simulators us- ing both physics based parameters and Neural Networks as universal function approximators. Our learning algorithm develops a mixed mode approach, mix- ing forward and reverse mode automatic differentiation with forward and adjoint based sensitivity analyses to efficiently perform gradient based optimization. We show that our physics informed learning method is capable of: (a) solving in- verse problems over the physically interpretable parameter space, as well as over the space of Neural Network parameters: (b) learning Lagrangian statistics of tur- bulence; (c) combining Lagrangian trajectory based, probabilistic, and Eulerian field based loss functions; and (d) extrapolating beyond training sets into more complex regimes of interest. Furthermore, our hierarchy of models gradually in- troduces more physical structure, which we show improves interpretability, gener- alizability (over larger ranges of time scales and Reynolds numbers), preservation of physical symmetries, and requires less training data.
INTRODUCTION
Why turbulence? Why Lagrangian? Understanding turbulent flows is crucial for many engineer- ing and scientific fields, and remains a great unresolved challenge of classical physics (Frisch) 1995). Many computational fluid dynamics (CFD) approaches (Harlow) (004) are being explored in this direction by approximating the Navier-Stokes (NS) equations which are well established in the field as explaining the "ground truth". One such CFD approach the mesh-free Smoothed Particle
Hydrodynamics (SPH) method (Gingold & Monaghan 1977| [Mnaghan 1992 2012). which has been widely applied to weakly- and strongly compressible turbulence in astrophysics and many engineering applications (Shadloo et al., (016)) It is one of very few approaches based on La- grangian construction: fluid quantities follow the flow using particles as opposed to the Eulerian approach which computes flow quantities at fixed locations by using computational mesh. This mesh-free, Lagrangian approximation of NS is appealing because it naturally unmasks correlations at the resolved scale from sweeping by larger scale eddies (Kraichnan 1964 1965). Developing approximation-opting SPH models, and then simulators, for turbulent flows is an ongoing area of research (Lind et al. 2020), to which this paper contributes. Specifically, we show how modern tools of machine learning and applied mathematics, such as deep neural networks (NNs), automatic differentiation (AD), and sensitivity analysis (SA), can be utilized and developed design better SPH models.
Physics Informed Machine Learning: What and Why? Numerical simulators, including classi- cal CFD simulators, have recently been blended with modern machine learning tools (King et al. 2018 Schenck & Fox 2018,, Mohan Gaitonde 2018, Mohan et al., [020a]] Maulik et al.) 2020 Ummenhofer et al., 2020 Mohan et al.) [200b]) Tian et al. (021) out of which a promising field, coined Physics Informed Machine Learning (PIML), is emerging (and being re-discovered (Lagaris et al. 1998)). As set, some eight years ago at the first LANL workshop with this name (CNLS at LANL. 2016, 2018, 2020). PIML was meant to pivot the mixed community of machine learning researchers one hand and scientists and engineers
the other, to discover physical phenom- ena/models from data. Today, PIML researchers are focused on incorporating known physical struc- tures (such conservation laws), physical hypothesis and physical equations (and/or their numerical

---

Under review as a conference paper at ICLR 2022
schemev/algorihhms) into modern machine learning tools. This development is exciting for physi- cists, as it integrates centuries of scientific discovery with modern machine learning - boosting new ideas and resurrecting interesting but forgotten hypotheses from the past. The overall goal of PIML is two-fold (King et al. 2018): learning new physics from data (experimental or computational) and improving/devellpingg machine learning algorithms by incorporating ideas from physics (especially in what concerns explainable, interpretable, and generalizable models).
Our Contribution: PIML for SPH. We develop learn-able hierarchy of parameterized "physics- explainable" SPH based fluid simulators that are trained and analyzed on SPH flow data. We train these models by mixing automatic differentiation (both forward and reverse mode) with forward and adjoint based sensitivity analyses. Utilizing these models and methods, we show that adding phys- ical structure improves interpretability, generalizability, and requires less training data compared to less informed models (such as Neural ODE). Additionally, we show the principal ability of these parameterized SPH simulators to solve inverse problems, and be fitted to the ground truth flow data using a combination of field based and statistical based loss functions. Furthermore, we show that the learned models are able to reconstruct the Lagrangian statistics of the flow, and can be used to learn unknown, missing, functions embedded within the SPH models by using Neural Networks as function approximators (Homik 1989).
Why PIML challenges are of interest to a broader Machine Learning Community? Machine Learning has made tremendous progress in recent years by developing a plethora of sound practical methods, algorithms and software which allow one to work with all kinds of data in application agnostic way. Today we can solve many challenging fitting and prediction problems of interpolation type, which were unthinkable even five years ago. However, there are still challenges that remain notabbly one related to extrapolation, and thus generalizability, into regimes where data is not sufficient, if available at all. One way to approach these challenges to substitute lack of data with cetiicatioo-specifffe modeling - in our case Lagrangian Modeling of Turbulent Flows. The emerging PIMI. community this research belongs to is focusing exactly on this. We are making progress, but also need help from a broader machine learning community on multiple issues, in particular these highlighted in this paper: How to embed physics-based into ML frameworks? How does embedding physics into ML schemes contribute (simplify or complicate) to roughness of the training landscape? How to balance Bayesian and deterministic approaches, e.g. selecting loss functions, in the intrinsically stochastic extrapolation problems, like turbulence? Finally, how to balance physics explanations (by tuning physically meaningful parameters, e.g. with sensitivity approaches) with predictive capabilities of machine learning (by training NN parameters, introduced for degrees of freedom which we do not need or do not know how to interpret)?
The overall outline of the remainder of the manuscript is a follows. First, in section 2, an overview of related studies is given that draws connections of this manuscript to the broader field of PIML. Next, in section 3] and Appendix A, formulation of SPH is presented, which provides the basis for our work on learning Lagrangian models of turbulence using parameterized and lear--alle SPH simulators. In section and Appenddx B. a detailed description of the learning algorithms, loss functions, and hierarchy of models are suggested, which provides the necessary tools to learn parameters of the SPH model. Then, in /ecctinn 5) and Appendix CJ we analyze this hierarchy of models and address how adding more of the known physical (SPH) structure affects the ability of the models to generalize to longer times and larger Reynolds numbers (Re). Finally, in sectionn 6 we draw conclusions on the progress made so far along with providing a discussion of future work.
Outline.
RELATED WORK
Early works of incorporating scientific domain knowledge (in the form of differential equations) within deep learning algorithms dates back to the 1990s (Lagaris et al.) 1998). However, in the context of modern deep learning, interest in this area has been revived (Raissi et al. 2017 Chen et al. 2019 Rackauckas et al. 2020 Ladický et al. 2015), due in part to the increased computational power afforded by parallelism across both CPUs and GPUs, along with notable achievements across disciplines (such scientific applications (Zhai et al., 2020), cota-comne algorithms (Wang jet al. 2016), computer vision (Serre) 2019), natural language processing (Li Deng 2018). etc. The main computational utilities and strategies of physics informed learning includes using physical
2

---

Under review as a conference paper at ICLR 2022
structure, NNs function approximators (Hornik 1989), automatic differentiation (Bücker 2006 and optimization tools to minimize an objective (loss) function.
There are several approaches that can be taken develop PIML algorithm: (1) using soft con- straints by adding physical structure, or symmetrics into the loss function (as was done in (Raissi et al. 2017) where Neural Networks are used discover and solve parameterized PDEs from data), (2) enforce physical constraints directly into the NN architecture (as (Mohan et al. 2020a) enforced the incompressibliitty constraint into Convolutional NN), (3) utilize NN's as function approxima- tors along with physical parameters/structur embedded within differentiable numerical simulators (as was done in (Chen et al., 2019; Rackauckas al., 2020 Tian et al.. (02T)), (4) mixing differ- entiable programming with SA along with physical structure and NNs as was done by (Ma et al.) 2021 Chen et al., 2019) (which from technical perspective is most related to our work). The broader focus of this work is on (4) as an approach to build parameterized Lagrangian models as candidates for SPH models of turbulence. However, this manuscript, we combine these scientific disciplines to not only provide tool for approaching the discovery of optimal Lagrangian models for turbulence, but also to explore the effects of adding physical structure into ML algorithms on extrapolation and generalizability.
In recent years, there have been several works that have integrated ML and DL techniques for La- grangian flows. The pioneering work by (Ladický et al. 2015) used SPH related models and regres- sion forests along with physics informed feature vectors, demonstrating that flow representations can be learned with data-driven techniques. In (de Anda-Suárez et al. 2018), evolutionary algorithms are applied for optimization of parameters in SPH flows. Further works began utilizing differentiable programming techniques, such as (Schenck Fox 2018) for robotic control for pouring applica- tions and in (Ummenhofer et al. 2020) where a continuous convolutional neural network operation developed on unordered particle data for learning and simulating SPH. (Tian al., 202T learns reduced models describing the Lagrangian dynamics of the velocity gradient tensor from Direct Nu-
merical Simulation data. Other works have used experimental data from Particle Image Velocimetry
and NNs to track particles embedding in flow (Rabault et al. 2017
Lee et al. 2017:
Cal et al.
2019
Stulov & Chertkov 202I).
SMOOTHED PARTICLE HYDRODYNAMICS
One of the most prominent particle-based Lagrangian methods for obtaining approximate numerical solutions of the equations of fluid dynamics is Smoothed Particle Hydrodynamics (SPH) (Mon- aghan 2005). Originally introduced independently by (Lucy) 1977) and (Gingold Monaghan| 1977] for astrophysical flows, however, over the following decades, SPH has found much wider range of applications including computer graphics, free-surface flows, fluid-structure interaction, bioengineering. compressible flows, galaxies" formation and collapse, high velocity impacts, ge- ological flows, magnetohydroddnamiic. and turbulence (Lind et al.) [0200] (Shadloo et al.) (0166). Below, we give brief formulation of SPH and in subsection 4.1 we hard code SPH structure into hierarchy of ppysiies-xpplanable" models.
3.1 APPROXIMATION OF EQUATIONS OF MOTION
Essentially, SPH is discrete approximation to a continuous flow field by using a series of discrete particles interpolation points (using integral interpolation with smoothing kernel W). Using the SPH formalism (see AAppendix A for more details), the partial differential equations (PDEs) of fluid dynamics be approximated by system of ordinary differential equations (ODEs) for each particle (indexed by i)
dr: dt U. dt
dv:
N
P3 P5
Σ
+ 1T;
V.W + Sexts Vi
(1)
See Appendix A for more detailed derivation, as well as description of how particle density and pressure P. are computed as in Monaghan (2012). this manuscript, deterministic and stochastic external forcing fext explored (see Appendix A for more details), which provides the energy injection mechanism. We also utilize the popular formulation an artificial viscosity II which approximates in aggregate contributions from the bulk and shear viscosity (ax), Vueman-Richtmen
3

---

Under review as a conference paper at ICLR 2022
viscosity for handling shocks (3) Monaghan (2012); Morris et al. 1997), as well the effective, eddy viscosity effect of turbulent advection from the under-resolver scales, i.e. scales smaller than the mean-particle distance.
3.2 EXAMPLE FLOW: TRAINING SET
Fig. ① shows consecutive snapshots of an exemplary multi-particle SPH flow in three dimensional space, where the coloration added for visualization purposes. We use a standard set of parame- ters for reakly-compressibe flows (see Cossins (2010))
1.0 (bulk-shear viscosity). 3 2a (Nueman- WWueman-Richtmyer viscosity) C = (speed of sound), ? = 7.0 with energy injection rate 0 0.5, and deterministic external porcing.(Appendie A). The inverse problem we pose consists in, given a sequence of snapshots, to reconstruct as best as we can the underlying (SPH) model used to generate the data.
(a) to
(b) 4
(c) l2
(d) ts
Figure SPH particles advancing in time used in training data, where coloration is added for visualization purposes and 4++1 4 R 100
MIXED MODE PHYSICS INFORMED MACHINE LEARNING METHODS
This section gives brief overview of our mixed mode physics informed learning algorithm (a first order gradient based optimization approach). This includes the hierarchy of parameterized SPH informed models, methods for computing the gradients, and formulating the loss functions. These methods are implemented using the open source software Julia (Bezanson et al. (017), which is a fast and flexible dynamic language, appropriate for high performance scientific and numerical computing, along with wide range of data science and machine learning packages such Flux, jl IInes et al., 22018), and AD packages Zygote.jl (Innes] 2018), ForwardDiff.jl (Revels et al., 2016). For more detailed treatment on learning algorithms (including derivations), and loss functions see Appendix B
4.1 HIERARCHY OF MODELS
In this work, we hard code SPH structure into hierarchy of parameterized models that includes physics based and Neural Network based parameters. Multilayer Perceptrons (MLPs) are used as universal function approximators (Hornik] 1989) embedded within an ODE structure describing the Lagrangian flows. It was found through hyper-parameter tuning that 2 hidden layers were sufficient for each model using a NN.
NODE: In this least informed model that explore, the entire acceleration term, that is the entire right hand side of Eq. (1) with the exception of fexts is approximated by an MLP, which is most related to the work done by (Chen et al. 2019). We make additional modification by considering the interaction of particles to be within a local cloud (using the cell linked list algorithm (Domínguez) al. 201 Note here that no pairwise interaction between particles is assumed. We assume that velocities, (t), and coordinates, (t), of N particles evolve in time according to
dr; du: dt
U. di
NNo (r: T. U2 1. Vj Ir. 1r.-,,1
< 2h) fext -i= 1,
N
(2)
where NN, 2ddm R' + RR 2d (d 2, is the space dimension, m is the fixed number of particles that closest to the i-th particle in each cloud, and I is the height, or number of
4

---

Under review as a conference paper at ICLR 2022
nodes, of the hidden layer). Although this NN is approximating a function which is interpretable (acceleration), the individual parameters the NN are not.
NN summand: To include more of the physical structure, represented via sum the i-th particle neighborhood Eq. ① (which is embedded in generating the "synthetic" ground truth data), we approximate the summand term from SPH using NN. Now pairwise interaction between particles is assumed. Here, Eq. ① is modeled by the Lagrangian based ODE
dr: dt
du,
Vi+ dt Σ NNo (r: Ty. - -;)) + fext Vi =i= N
(3)
where NNO: : R2d + R' - R R4.
Rotationally Invariant NN: In this formulation, built on the top of the NN summand, we use neural network of the form NN I - + + R to approximate the pair-wise part of the acceleration term in Eq. ( where the rotational invariance is hard coded by construction (about simple rotationally invariant basis expansion using the difference vector (r; r2)), as follows
dr: = u, dt dt
N dv: = Σ NNo
P. P3 (r. rg) ..) (v. - 1,). -ryll (r: - ry) + fezt (4) p p
VP. NN: Next, we push the neural network deeper within the summand term and now explicitly include the II, i.e. artificial viscosity term, and use a NN approximate the pressure contribution (i.e PP term in SPH Eq. ① see Appendix for more details):
dr:
dv, Σ (NNO (r: -r) ry) VWy) + fext =1=1... ,N
Uis dt
(5)
dt
where NN g: : Rd R R
EoS NN: Embedding NN within an SPH simulator for learning the interpretable equation of state from flow data using a neural network (Prig) for approximating P(p) in Eq. (I:
dr:
dv; dt
(pi)(pp)
Pnno(ps))
Us
Σ
+ Ili VW + feet Vi -i= 1 N
(6)
dt
P
P
where Puno (p)
R +R
Fully Physics Informed: In this formulation, the entire physical structure, that is known to gen- erate the ground truth, is used, and the physically interpretable parameters ax, B, 7, e are learned
dr; = vis dt
dv,
P.(c,))
+ Py(e,)) -,y(c,,, VW + fext -i= 1, N. p
Σ
dt
p
We note that as more of the physical structure is added into the learning algorithm, the learned mod- els have more interpretability; the learned parameters associated with actual physical quantities.
4.2 MIXING SENSITIVITY ANALYSIS AND AUTOMATIC DIFFERENTIATION
Sensitivity analysis (SA) is a classical technique found in many applications, such as gradient- based optimization, optimal control, parameter identification, model diagnostics, (see (Donello et al.) 2020) and many historical references there in), which was also utilized recently to learn neural net- work parameters within ODEs (Chen al.) 2019 Rackauckas et al.) 2020 Ma al. (021). In the context of this work, we use SA to compute gradients of our hierarchy of parameterized models. We mix SA with Automatic Differentiation (AD) (see (Ma al., 2021 Bücker 2006) and references there in): forward mode and reverse mode AD is applied to derivatives within the SA algorithm, where the method chosen based on efficiency (depending on the dimension of the input and out- put space of the function being differentiated). The losses used in this work are defined in Eq. (4.4), but in general, are, respectively, the vector of the combined vector of the particles" coordinates and velocities and the vector of parameters introduced in the next subsection.
consider loss functions of the form, L(X, fe ((X 0, t)dt, where X and
5

---

Under review as a conference paper at ICLR 2022
4.2.1 FORWARD AND ADIOINT BASED METHODS
Let us, first, introduce some useful common notations for our SPH informed models, X; (ri, u,)", X = {X,]i I,,.N), where, (Ti: Vs. 21), and, = (uz, 2g W/), are the position and velocity of particle i respectively. Also, = , OP|2 where p is the number of model parameters. Now, the SPH discretizations defined in Eq. 1 can be stated in the ODE form
Vi: ((
(v,, E().
(8)
where F; the right hand side of Eq. I Forward and Adjoint based Sensitivity Analyses (FSA, ASA), analogous to forward and reverse mode AD respectively, are used to compute the gradient of the loss function (see Section 4.4 0,l(XX 0) f 0x (xx 0, 1)dd X ((, + V 0, l)dl, (for derivations of FSA and ASA see Appendix B FSA computes the sensitivities, S. = d by simultaneously integrating a system of ODEs:
Vi, Va = 1, .............................. .P: dS, (aF 0)/XX) + OF; 00/00
(9)
which is known to be more efficient when the number of parameters, p. is not very large, S (100)). Once S is known, the gradient is computed directly from the formula for ôg L above. The ASA avoids needing to compute de X by instead numerically solving system of equations for the adjoint equation backwards in time, according to subsubsection B.0.2 Once is found, the gradient of the loss function is found according to, 0gL = .- Σ. A ( dt.. It is well known that the ASA more efficient when 2 (1100) (see (Ma al., 2021), D), because solving the adjoint equation is independent of p but requires more memory to store forward the solution then to integrate A backwards in time.
4.3 MIXED MODE AD
In both FSA and ASA described above, the gradient of F; with respect to the parameters, OF, (X)),, (T))0))00 0)/00, and the Jacobian matrix, {aF, 0)/00)/iiii,... need be computed. In this manuscript, we accomplish this with a mixed mode approach, i.e. mixing forward and re- verse mode AD, where the choice is based on efficiency. Depending on the model used, many of the functions to be differentiated have varying input and output dimensions. For example, when 2d if p >> 2d, then reverse computing aF, (x (T), 0)/08*, with AD where, (())) R"
mode AD is more efficient than forward mode (Bücker 2006).
4.4 Loss FUNCTIONS
In this section, we construct three different loss functions: trajectory based (Lagrangian), field based Eulerian, and Lagrangian statistics based, described the following three subscctions. Since our overall goal involves learning SPH models for turbulence applications, it is the underlying statistical features and large scale field structures we want our models to learn and generalize with. This is discussed further in section 5 where compare the hierarchy of models that constructed in subsection 4.1
4.4.1 TRAJECTORY BASED Loss FUNCTION
A simple loss function to consider is the Mean Squared Error (MSF) of the difference in the La- grangian particles positions and velocities, they evolve in time, L,r = MSE(X, (()) |X (0))12// where X and X are the particle states the ground truth and the predicted, respectively. Minimizing this loss function will result in discovering optimal parameters such that the predicted trajectories gives the best possible match (within the model) for each of the particles. Notice that this "perfect" matching of the multi-particle state is not really appropriate for tracking temporal correlations in turbulence on the time scales longer than the turn-over time of the resolved eddy (estimated, roughly, the time needed for pair of initially neighboring particles to separate on a distance comparable to their initial separation). This is due to the fact that turbulence intrin- sically chaotic, therefore resulting in strong sensitivity of the state (Lagrangian particle positions and velocities) on its initial conditions. Therefore, we expect the trajectory based loss method to over-fit fit when the observation (tracking) time is sufficiently long.

---

Under review as a conference paper at ICLR 2022
4.4.2 FIELD BASED LOSS FUNCTION
The field based loss function tries to minimize the difference between the large scale structures found in the Fulerian velocity fields, [0(0)= = MSE(V), V = I //I//N where V (mg/
W -ryll, h) uses the same SPH smoothing approximation to interpolate the partiele velocity onto a predefined mesh r (with Ns grid points). Let us remind that SPH is, by itself, an approximation for the velocity field, therefore providing a strong additional motivation for using the field based loss function.
4.4.3 LAGRANGIAN STATISTICS BASED LOss FUNCTION
In order to approach the statistical nature of turbulent flows, one can use well established statis- tical tools/objects, such single particle statisties. In this direction, consider the time integrated Kullback-Leblerr divergence (KL) as a loss function
Ege((), ( Ggell, (i) = S. J- Gtelt, (t)(t), log (pread (1, 2pr (0, 1), x)) dxdtt
where (t), and (0, represent single particle statistical objects over time of the ground truth and predicted data, respectively. For example, we can use the velocity increment, z x'()) = (ōu,, Öw;,),;, where (t, (t) (0) and z ranges over all particles. Here G(t, x(t). x) is a continuous probability distribution (in 2)) constructed from data x(4) using Kernel Density Es- timation (KDE), to obtain smooth and differentiable distributions from data, as discussed in (Chen) 2017). that is G(T, x) (Nh) -1 Σ K ((z; - =)/h), where K is the smoothing kernel (cho- sen to be the normalized Gaussian in this work). (See Appendix [B. for details.) In 3D experiments (discussed in the next section), we use Lu + Ly as a combination of statistical and field based loss functions, with the expectation that the gradient descent (in the process of training) will drive different parameters, responsible for large and small structures, to their optimal values unison.
5 RESULTS: EVALUATING AND TESTING THE HIERARCHY OF MODELS
In this section, we show the ability of the parameterized SPH simulators to; solve inverse problems, fit underlying flow data using a combination of field based and statistical based loss functions, learn unknown functions embedded within SPH based models, generalize to flows not seen in training and compare their ability to conserve linear and angular momentum. In subsection [.] we illustrate the ability of different models in the hierarchy to learn the parameters to approximate the (underlying) distribution of velocity increments. Then, in section 5.2 we focus on comparing the generalizability of different models in the hierarchy. We show that the more physical structure is hard-coded into the model, the better is the generalizability ( (i... extrapolation over larger time scales and Reynolds numbers) and the better it conserves linear and angular momentum (preserves transnational and rotational symmetries respectively). We also notice that another important advantage of the physics informed models is that they require less training data (as seen in Figure 9).
5.1 SOLVING INVERSE PROBLEMS
Each parameterized SPH model within the hierarchy (see subsectiinn 4.1) trained under equivalent conditions: (a) on the same SPH samples (see subsection 3.2) of fixed temporal duration (which we choose be equal to the time scale required for pair neighboring particles to separate by the distance which is in average factor O(1) larger than the pair's 's initial separation and henceforth de- noted fx): (b) with the same loss function; and (c) with a deterministic forcing foxt (subsection A.3) with constant rate of energy injection.
In Figure 2 and Figure 6) we see that the physics informed parameterized SPH simulator is learn- able; the physical parameters a, B, c, 7 are learned over the physically interpretable parameter space. Figure 3 jillustrates the ability of the method applied to the EoSNN model to learn (approximate) physically interpretable functions (P(p)) using NNs embedded within SPH model. FFiuure 4 shows that each model is capable of learning the parameters so that the underlying velocity increment distribution (a single particle I Lagrangian statistic) is approximated, i.e learning Lagrangian statistics.
7

---

Under review as a conference paper at ICLR 2022
These figures, show that mixture of loss functions described in subsubsection 4.2.1 the SPH informed models subsection 4.1
subsection 4.4 applying the methods
are learn-able.
WC.SRL
WC_SPH_ AV: N 4096, time
18
12
Figure 2: Solving an inverse problem for the fully physics informed model on 3D SPH flow with deterministic external forcing on 4096 particles over 4 SPH parameters. The solid lines show the SPH model parameters (initially chosen to be uniformly distributed about (0, ()) converging to the dashed lines representing the ground truth parameters. Here the field based loss function (see (sub- (section 4.4) is used and is converging up to the order of machine precision (as Float64 and the MSE of field difference was used). See Figure and Figure for similar results.
Learning 505 WCSPH N 4096, T
PMLL SPHLAV Wrandn N 4096, time steps
12"
200
3000
tsraloons
Figure 3: Using a Neural Network to approximate the equation of state (see subsection 4.1 using the field based loss function. We see that P(p) is well approximated, where Pgt is the ground truth EoS and Pun is the neural network approximation. We note that the underlying ground truth data has density variations within 1% of mean density, so the NN sufficiently approximates the EoS over a larger domain of densities as seen in training (see Figure for another example).
1A
12
KALAA
Figure 4: The velocity increment distribution (in u, i. G&n. learned over each model in the hierarchy subsection 4.1 The dashed crossed green line is the predicted velocity increment distribution the specific iteration Gg. and the solid purple line is the ground truth velocity increment distribution. Each model is learned using the Lu + Ly loss and the FSA method. From left (NODE) to right (Fully informed PhysIn fnf) more physical structure added. See Figure 5) for further comparisons of these learned models
5.2 GENERALIZABILITY (EXTRAPOLATION CAPABILITY)
When the training is complete (i.e. when the loss function reaches its minimum) we validate ex- trapolation capability of the models test data sets which longer in duration derived from the setting corresponding to stronger turbulence (which we control by increasing intensity of the injection term, fext; while keeping the integral, i.e. energy injection scale, constant).
8

---

Under review as a conference paper at ICLR 2022
Spacnday
187
301
50
355
0% Dff
53%
(a) t Generalizability
(b) Re Generalizability
(c) Rotational error
(d) Translational error
Figure 5: (a) Box plot using 15 different generalization errors over physically relevant time scales ranging from the shortest, l (which is the time it takes for a pair of initially neighboring particles to double their relative separation), up to the longest, leddy (which the time takes for the largest eddy of the flow to make complete revolution). (b) Box plot over 5 different generalization emrors over different Reynolds numbers, ranging in values up to 3 times larger than training data where the time scale is fixed at xx. (c) Error in rotational and (d) translational symmetries in the learned models (i.e measuring conservation of angular and linear momentum of learned models respectively).
We observe in Figure 5| structure added to the model is, the lower the generalization error becomes when we extrapolate to flows with larger Re numbers (and over different time scales). Further experiments were conducted on two dimensional flows (see Figure 10, shown in the Appendix; here, we observe similar trends in generalization). Furthermore, Figure 5) shows that as more physical structure is added to the model, the error in rotational and translational invariance decreases (and therefore conservation of angular and linear momentum is preserved, even though it was not enforced directly). Finally, we additional qualitative generalizability results in Figure 9 (shown in the Appendix) which allows us to conclude that data (samples of particle trajectories) of much shorter (temporal) duration are required to train models which are more physies informed (where a short time duration, (x, is used in training, and we see the converged models extrapolating to longer time scales not seen in training). Figure also shows, qualitatively, that the large scales structures of the flows are well preserved as we extrapolate for longer times than used in training.
showing results for trained models using Ly + Lit, that the richer (physical)
6 CONCLUSIONS
Combining modern tools in CFD, machine learning, deep learning, automatic differentiation, and classical sensitivity analysis, we have developed a learn-able hierarchy of parameterized "physics- explainable" Lagrangian fluid simulators and showed that adding physical structure improves inter- pretability, generalizability (over different time scales and Re numbers), and requires less training data compared to less informed models (such as Neural ODE). Additionally, we see that as more of the physical structure is hard codded into the SPH based models, the better is its ability to conserve linear and angular momentum (translational and rotational symmetries respectively) which is known to be the case with the underlying ground truth flow.
Furthermore, these parameterized SPH simulators can be used to solve inverse problems (in the physically interpretable parameter space as well as in the NN parameter space), fit underlying flow data using combination of field based and statistical based loss functions, and be used learn unknown functions embedded within SPH based models. We showed that each model is capable of learning the parameters so that the underlying velocity increment distribution (a single particle Lagrangian statistics) is approximated.
In the future, we plan to improve beyond what is accomplished in this manuscript, that is solve efficiently inverse Lagrangian problems with SPH. Specifically, we aim to build Reduced-Orie La- grangian Models of Turbulence based on the SPH approach which provide the best to turbulence data coming from more realistic simulations and/or experiments. It is of special interest to push the PIML based Lagrangian SPH approach towards discovering reduced turbulence models which will extrapolate into difficult regimes, e.g. corresponding to regimes of high Re number and stronger degree of compressibility, where the data are limited.

---

3-17778--3319770000000 Under review as a conference paper at ICLR 2022
7 REPRODUCIBILITY STATEMENT
Supplementary material is provided that includes the Julia source code of our parameterized SPH simulators, the gradient based learning algorithm, and sensitivities for each model in our hierarchy. The general steps to take in order to reproduce these results are as follows: generate SPH ground truth flows using the parameters described in section 3 (or use the data-sets provided), then use this data for learning the hierarchy of models found in the main.jl file (under the 3d-phys-semi-in inf direc- tory) which requires a selection of the model, loss function, sensitivity method, number of iterations, and height (where height is found through tuning to have value of 5, other than for the EoS model which performs better with a height greater than 8. Note, with these heights, most of the models containing NNs have total number of parameters usually fastest (ASA is slightly faster for the NODE model, as it involves more parameters than other models due to non-pairwise interaction of particles). However, for consistency, we have selected the FSA method when comparing all models with this relatively small number of parameters. The saved models are provided in the output data files which are used in processooss.igg Furthermore, in the final version of this manuscript a link to the github repository will be included.
the order p (110)); in this case FSA is
REFERENCES
Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral Shah. Julia: fresh approach to nu- merical computing. SIAM review, 59(1):65-95) 2017. URL
1137/
141000671
Andrew Bradley. PDE Constrained Optimization and the Adjoint Method. cs.stanford.edd 2019.
URL //ttts://cs....... wtw.tord.eddu wonod//dttootntttttitalllll pdf
Martin Bücker. Automatic Differentiation Applications, Theory and Implementations. Springer, 2006.
Shengze Cai, Jiaming Liang, Qi Gao, Chao Xu, and Runjie Wei. Particle image velocimetry based on deep learning motion estimator. IEEE Transactions Instrumentation and Measurement, 2019.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary dif-
ferential equations. arXiv: /806.0736066 2019.
Yen-Chi Chen. tutorial on kernel density estimation and recent advances.
Biostatistics & URL https:
Epidemiology, 1(1):161-187, 2017, doi: 10. 10.1000/2470956000117..333
7/doi. .rg.rrrgggo/ 10. 1010000000300000 2017 7 1396742
CNLS LANL. Physies Informed Machine Learning (PILL) workshop.
http:
/ //ww. event t. com httpc///eeeeeeeeaaaaaaaaaaaaaaee
event raa.maaryryr89grr http://www...2333333666666666660
2016, 2018, 2020.
Peter J. Cossins. Smoothed particle hydrodynamics. arXiv:1007.1245,000 2010.
Juan de Anda-Suárez, Martín Carpio, Solai Jeyakumar, Héctor Soberanes, Juan Mosiño, and Laura Cruz Reyes. Optimization of the Parameters of Smoothed Particle Hydrodynamics Method, Using Evolutionary Algorithms, pp. 153-167. Springer, 01 2018. ISBN 978-3-319-71007-5, doi: 10.
J. M. Domínguez, A. J. C. Crespo, M. Gömez-Gesteira, and J. C. Marongiu. Neighbour lists in smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids, 67(12):2026-2044 2042, 2011. doi: https://doi.or//0..ooo 24 URL https:// enlinelibrary wat.ary.willyy.. com/ doi/ abs/100.002/Eldd.
Michael Donello, Mark Carpenter, and Hessam Babaee. Computing sensitivities in evolutionary systems: A real-time reduced order modcling strategy. arXiv:2012.14028, 2020.
U. Frisch. Turbulence: The Legacy of AN Kolmogorov. Cambridge University Press, 1995.
10

---

Under review as a conference paper at ICLR 2022
David A. Fulk and Dennis W. Quinn. An analysis of 1-d smoothed particle hydrodynamics ker- nels. Journal of Computational Physics. 126(1):165-180) 1996. ISSN 0021-9991, doi: https://
1996.0188. URL /rttps://ae www.sciinne.iettttt..mmmmm wttp://iiiettt.cccccoceoooata art.cle/p1i//0000 9991 9690128X
R. A. Gingold and J. Monaghan. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 181:375-389, Nov 1977.
Francis H. Harlow. Fluid dynamics in Group T-3 Los Alamos National Laboratory: (LA-
UR-03-3852). Journal of Computational Physics, 195(2)4414-33,, 2004. ISSN 0021-9991. doi: htips://doi.co com/ http:///ceoooaecatcccee/////./..
URL htttp::/wwww.cienneee
et al Hornik, Kurt. Multilayer feedforward networks are universal approximators. Neural Nenworks, (15):359-366, 1989.
Michael Innes. Don't unroll adjoint: Differentiating ssa-form programs. CoRR, abs/1810.07951, 2018. URLAttp::/arxxtt wwtp:/iwwwr.aa.aa..a./aa/ap.o.aa /wt...111.... 08901
Michael Innes, Elliot Saba, Keno Fischer, Dhairya Gandhi, Marco Concetto Rudilosso, Neethu Mariya Joy, Tejan Karmali, Avik Pal, and Viral Shah. Fashionable modelling with flux. CoRR, abs/1811.01457 2018. URL
Ryan King, Oliver Hennigh, Arvind Mohan, and Michael Chertkov. From Deep to Physics-Informed Learning of Turbulence: Diagnostics. arXivv1810.07785, 2018.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv:1412.6980, 2017.
Robert H. Kraichnan. Kolmogorov's hypotheses and eulerian turbulence theory. The Physics of Fluids, 7(11):1233-1734, 1964. doi: 10. 174677... URL httthtt/p//:ppp / sciiaattioo......
org/doi/abor 10. 1et.1063111.0 2746572
Robert H. Kraichnan. Lagrangian-histee history closure approximation for turbulence. The Physics of Fluids, 5(4):575-598, 1965. doi: 10. 176217 URL
org/dol/abs/11 1063/1..0......
L'ubor Ladický, SoHyeon Jeong, Barbara Solenthaler, Mare Pollefeys, and Markus Gross. Data- Driven Fluid Simulations Using Regression Forests. ACM Trans. Graph. 34(6), October 2015.
ISSN 0730-0301. doi: 10. 19.1145/281699528888222
URL
2816795.0 2818129
I.E. Lagaris, A. Likas, and D.I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 95),987-1000, 1998. ISSN 1045- 9227. doi: 10.1109/72.712178. URL httpt:::://xoooo//////tt
doi. org. wtt.://///......//..... 1109/72. 712178
Yong Lee, Hua Yang, and Zhouping Yin. Piv-dcnn: cascaded deep convolutional neural networks for particle image velocimetry. Experiments Fluids, 58 58(12):171, 2017.
Yang Liu Li Deng. deep learning in natural language processing. Springer, 2018.
Steven J. Lind, Benedict D. Rogers, and Peter K. Stansby. Review of smoothed particle hy- drodynamics: towards converged lagrangian flow modelling. Proceedings of the Royal So- ciety A: Mathematical, Physical and Engineering Sciences, 476/2241).0190880,, 2020. doi: 10. 10.10988/9pa.2019.088) URL bttpes///bopsoociitppppaaooooo.p
10 10..199/rsss hwt.9////pp.... 2019.....1..880
L. B. Lucy. numerical approach the testing of the fission hypothesis. Astronomical Journal, 82: 12:1013-1024, December 1977. doi: 10. 10.11066/12211
Yingbo Ma. Vaibhav Dixit, Mike Innes, Xingjian Guo, and Christopher Rackauckas. A compari- son of automatic differentiation and continuous sensitivity analysis for derivatives of differential equation solutions. /rXiv:18/2.0/800 2021.

---

Under review as a conference paper at ICLR 2022
Romit Maulik, Arvind Mohan, Bethany Lusch, Sandeep Madireddy, Prasanna Balaprakash, and Daniel Livescu. Time-series learning of space dynamics for reduced- order model closure. Physica D: Nonlinear Phenomena, 405:132368, Apr 2020. ISSN 0167-2789. doi: 10.1016/j.
physd 132368. URL http://tt//// /dbp::///xx. doi. wrg.orgg//... /0. 101 physd. 2020. 132368
Arvind T. Mohan and Datta V. Gaitonde. A deep learning based approach to reduced order modeling for turbulent flow control using istm neural networks. https//furkiioorrgaa///11111///1 2018.
Arvind T. Mohan, Nicholas Lubbers, Daniel Livescu, and Michael Chertkov. Embedding hard phys- ical constraints in neural network coarse-gras -ratniing of 3d turbulence. arXiv:2002.0002/, 2020a.
Arvind T. Mohan, Dima Tretiak, Misha Chertkov, and Daniel Livescu. Spatio-temporal deep learn- ing models of 3d turbulence with physics informed diagnostics. Journal of Turbulence, 21(9- 10):484 524, 2020b. doi: 4665480206 202...2233000 1080 1080/146888 4652448 2020. 832230
J. Monaghan. Smoothed particle hydrodynamics. Annual Review of Astronomy antrooy 30(1):543-574, 1992. doi: 10.1 10. 10.1146/anmurer/ aa 30 00.00999922.. 0001920002251
http:/rrww.......00992000055555. URL/attps://doi.org
Monaghan. Smoothed particle hydrodynamics. Reports on Progress Physics, 68(8):1703 1759, Aug 2005.
J.J. Monaghan. Sph and riemann solvers. Journal of Computational Physics, 136(2):28-3077, 1997. ISSN 0021-9991. doi: http:///o..rr//11106iighhh.....9 sciencedirect com/ htt.://cww...c..cooooocoooooaota
URL httt:ttttt///tt
J.J. Monaghan. Smoothed particle hydrodynamics and its diverse applications. Annual Review
of Fluid Mechanics, (4(1):323-333)) 34,,3333,,,, 2012. doi: 10.1 20.1146/amurev--fuuie 1201100-0220000001222200 URL
https://doi. www.o//......./.../e......
fluid 110110
J.J Monaghan and R.A Gingold. Shock simulation by the particle method sph. Journal
of Computational Physics, 52(2):374-383 1983. ISSN 0021-9991. doi: htps://doi.org!
001101600211-9999188399003333330 URL httpp:://////t www.ecce.e.e....e.rrett......... wwt.rreeetcccce.ee.e comm/ccience/ arricee/s11/002199918900
Joseph Morris, Patrick Fox, and Yi Zhu. Modeling low reynolds number incompressible flows using sph. Journal of Computational Physics, 136(1):114-266, 1997. ISSN 0021-9991. doi: https://doi.org/I 1006///pp............. 1997.5776. URL htttps://ww...ttt. /twp:///wwwww....oo http///wwwwieeeeeeeeettttttt com/
M. Petersen and D. Livescu. Foreing for statistically stationary compressible isotropic turbulence. Physics of Fluids, 22111):11101 2010.
Jean Rabault, Jostein Kolaas, and Atle Jensen. Performing particle image velocimetry using artifi- cial neural networks: a proof-of-concept. Measurement Science and Technology, 28(12):12530) 2017.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scien- tific machine learning.
04385 2020.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations. arXiv:1711.10566, 2017.
Revels, M. Lubin, and T. Papamarkou. Forward-mode automatic differentiation in Julia. arXiv:1607.07892 [cs.MS], 2016.
Connor Schenck and Dieter Fox. Spnets: Differentiable fluid dynamics for deep neural networks, In Aude Billard, Anca Dragan, Jan Peters, and Jun Morimoto (eds.), Proceedings of The 2nd Confer- on Robot Learning, volume 87 of Proceedings of Machine Learning Research, pp. 315-335. PMLR. 29-31 Oct 2018. URL ht html
mir. press/v8 /ttsss///8/sschncckkaa. sch.ncceec..aacc
12

---

Under review as a conference paper at ICLR 2022
Thomas Serre. Deep learning: The good, the bad, and the ugly. Annual Review of Vision Science, (1):399 426, 2019. doi: 10. /0.1146//mmurevvvvssioon 091718 8-014951. wwtp://aaaaaeeeiiooooo---o------
org/
PMID: 31394043.
M.S. Shadloo, G. Oger, and D. Le Touzé. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges. Comput- & Fluids, 136:11-34, 2016. ISSN 0045-7930.
doi: https://bb..rr//1100066....opppp
2016.05.029 URL https://wwes http://wwwwwccccmm
50C4579301 6301773
Nikolay Stulov and Michael Chertkov. Neural particle image velocimetry. arXi::2101.11950 2021.
Yifeng Tian, Daniel Livescu, and Michael Chertkov. Physics-informed machine learning of the lagrangian dynamics of velocity gradient tensor. Phys. Rev. Fluids, 6:094607, Sep 2021. doi: 10.1 20.110//PyysReeFuudd... 6.094677. URL http:s///..nka.s..rrg////// PhysRevF wwwds 6. 094607
Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian Fluid Simu- lation with Continuous Convolutions. In International Conference Learning Representations,
2020.
Yasi Wang, Hongxun Yao, and Sicheng Zhao. Auto-encoder based dimensionality reduction. Neu- recomputing, 184:232-242 2016. ISSN 0925-2312. doi: https://doi.org/10...00 neucom. 2015.08. 104. URL https://wwe/ http://www..iccecec/ee encedi rect. com soon/scientng htt///ciinne/aarreee clt.tii//ccp////// $092523121501766 RoLoD: Robust Local Descriptors for Computer Vision 2014.
www.sccienediii
Xiaoming Zhai, Yue Yin, James W. Pellegrino, Kevin C. Haudek, and Lehong Shi. Applying ma- chine learning in science assessment: a systematic review. Studies in Science Education, 56 (1):111-151, 2020. doi: 10.108000305726..2200..777 1080/03057267. 2t1p:/3305227.2222.112222
URL https: //oo///o//... lorado..orrooo //tp:///o/o....../////
A APPENDIX: SPH (IN MORE DETAILS)
A.1 BASIC FORMULATION
Essentially, SPH is discrete approximation to a continuous flow field by using a series of discrete particles. Lets first start with the identity
A(r) = A(r))((-- A
(10)
where is any scalar or tensor field. Using the smoothing kernel W (for interpolation onto smooth "blobs" of fluid) and after a Taylor expansion it can be shown (according symmetry of smoothing kernel [Cossins (2010) that
A(r) - -/-
(11)
where Wi is constrained to behave similar to the delta function,
J W hydr =ddr111 = lim W(r, h) = 5(r). h-+
The choice of smoothing kernels is important, and effects the consistency and accuracy of results Monaghan 2005). where bell-shaped, symmetric, monotonic kemnels are the most popular Fulk & the best smoothing kemnels to use. Commonly
Quinn (996). however there is still disagreement used are the B-spline smoothing kernels with a finite support (approximating a Gaussian kernel). The cubic smoothing kernel is used in this work has the following form:
1 (2 9)3 (1
0 0 q <
(())
(2
q)3
1 q 2
(12)
2 2
13

---

Under review as a conference paper at ICLR 2022
where, W(r h)l,h) h w(9) with q |r ryl/h. and
o(d)
[1/T if d = 3. 10 if 2 a normalizing constant to satisfy the integral constraint on W (see Cossins 2010); Fulk & Quinn (1966); Monaghan (2012) for more details). The finite support allows one to use neighhorhood list algorithms discussed below to utilize computational advantages
(only requires local cloud of interacting particles instead of all the particles in the computational domain that would be required with a Gaussian kernel).
SPH be formulated through approximating integral interpolants of any scalar or tensor field A by a series of discrete particles
(A(r)) = /
Σ
A(r,) m, p(r,)
(1)
(13)
(i.e a convolution of A with W) where dr' denotes a volume element and W(q. h) is the smoothing kernel. Each particle represents continuous "blob" of fluid and carries the fluid quantities in the Lagrangian frame (such as pressure P;, density Ps, velocity U; etc.)
The convenience of this method becomes apparent when the differential operators are approximated, i.e.
V, (A(r)) dr /, (1r-)22,, s0
A(r Σ m V W (11l)1)),)))))))
(14)
where we see that in this direct approach to approximate the gradient operator we only need to know the gradient of the smoothing kernel (which usually fixed beforehand). Multiple methods have been proposed and different methods are best suited for different problems. Similar approximations hold for taking the divergence or curl of vector field Gingold & Monaghan (1977). The most common approximation of the gradient operator is derived from the following identity,
v,A(r:) = (A
A(rr)) V
(A(Ti)
Vr
(15)
and is approximated with particles as
N
A
A
V,A Ps
mj
Vr Wag.
(16)
A.2 APPROXIMATION OF FLOW EQUATIONS
The above integral interpolant approximations using series of particles can be used discretize the equations of motion (as seen in Eq. (I) with more details found in GGigold & Monaghan (1977); Monaghan 2012). Each particle carries mass m, and velocity v, and other properties (such as pressure, density etc.). We can use Eq. (13) to estimate the density everywhere by
p(r:) = Σ my (m:
(17)
where although the summation is over all particles, because the smoothing kernel has compact sup- port the summation only needs to oceur over the smoothing radius (here 2h. as seen in Eq. (12) Another popular way approximate the density is through using the continuity equation and ap- proximating the divergence the velocity field in different ways|Monaghan (1992). In what follows we use the notation A; A(r:). Using the gradient approximation defined above (Eq. (16), the pressure gradient could be estimated by using
p V. Σ my(P, P.)V, W.j.
where Way -ryl, However, in this form the momentum equation d,v --V does not conserve linear and angular momentum Monaghan (992). To improve this, a symmetrization is
14

---

Under review as a conference paper at ICLR 2022
often done to the pressure gradient term by rewriting V.P = 0r (2) + _ vnp. This results in momentum equation for particle discretized as
du: dt
P P
Σ m4 3
V Wag
which produces a symmetric central force between pairs of particles and result linear and angular momentum are conserved Monaghan (1992). Next, including an artificial viscosity term and external foreing, the full set of ODEs approximating PDEs governing fluid motion is
dr; dt
=
Vi f1,2....NJ
(18)
dv:
P3
P: P
Σ j71
m
+ Hنi
rr Wi + fent Vi (I.
(19)
In this work, start by using the weakly compressible formulation by assuming barotropic fluid, where equation of state (EoS) given by
P(p) = 2:P0 7
( Pb
-1
-1
(20)
as Monaghan (2012), where Po is the initial reference density, and y = is used. In future work, we plan on including the energy equation to extend these methods for highly compressible applications.
There are many different forms of artificial viscosity that have been proposed Monaghan (997)) In this work, use the popular formulation of I1, that approximates the contribution from the bulk and shear viscosity along with an approximation of Nueman- -ichtmye viscosity for handling shocks Monaghan (2012); Morris al. (1997):
axeet + 3p4
Fty <
Pi
(21)
0.
otherwise
where è represents the speed sound and
hvi rgg
Maj
r.y 2 + ch2
Pg 0.5(pt(p: Pp),
This artificial viscosity term was constructed the standard way following Monaghan 1997); Mon- lghann & Gingold 1983): The linear term involving the speed of sound was based the viscosity of an ideal gas. This term scales linearly with the velocity divergence, is negative to enforce II, > O, and should be present only for convergent flows (V.j Tú < 0). The quadratic term including (vij T2 cial pressure roughly proportional to 2/ and approximates the von Neumann- Richtmme viscosity (and should also only be present for convergent flows). There are several advantages to this formu- lation of Ily; mainly it is Galilean and rotationally invariant, thus conserves total linear and angular momentum. more detailed derivation is found in Cossins (2010) (along with other formulations of artificial viscosities).
2 is used prevent penetration in high Mach number collisions by producing
artifi-
In practice the summation Eq. over all particles is carried out through a neighborhood list algorithm (such as the cell linked list algorithm with computational cost that scales as O(N) Domínguez et al. (2011)). We also note that Eq. (T) can also be derived from Euler-Lagrange equa- tions after defining a Lagrangian, see Cossins (2010) (for respective analysis of the inviscid case, when the artificial viscosity term is neglected), then an artificial viscosity I1, term can be incorpo- rated from using SPH discretizations, see Monaghan 2005) for details.
15

---

Under review as a conference paper at ICLR 2022
A.3 EXTERNAL FORCING
In order to approach a stationary homogeneous and isotropic turbulent flow, we use a combination of external forcing, both stochastic and deterministic. We first stochastic external forcing acting on large scales which forces 10 modes that are close to the origin in Fourier space. Formally,
fegt Sext
Sezt
Seal = -E== cos cos [ + JK + 2T(((
2.5
v(FE)² + (fE)
Se [1] = -Et 10
cos(JEa: + fluv + 2n((t, β = 0.25
1.5
where (fk: fk) are the frequency components that lie near to the origin (in frequency space), and C(t, k) is stochastic white noise term.
Once the flow achieves stationarity, deterministic forcing is used (for simplifying the learning algorithms described the next section), which is commonly used CFD literature, e.g. (Petersen & Liveseu 2010) for analyzing stationary homogeneous and isotropic turbulence. Then,
ke (v: B), ke 0.5
N u
fext [
+ +
the kinetic energy computed at each time step, 0 represents the rate of energy injected into the flow and is the mean velocity for each component. is fixed quantity which is set before hand, and the flow is integrated until stationarity is reached.
A.4 NUMERICAL ALGORITHM FOR FORWARD SOLVING SPH
We use a symplectic integrator for generating the "synthetic" ground truth data, and for making prediction steps required in our gradient based optimization described in Eq. (4 We apply the Verlet integration scheme (leap-frog) to solve the Initial Value Problem. Using the notation,
dr: = vi dt
dv: dt
X = ((ri, vi {1,... N)}, p = (pilii (1....,N)) N}}.
F.(p. X),
we proceed according to the following algorithm
1: Compute p using Eq.
13
2: Compute F. (p", X²) using Eq. 1)
3: v;
k+1
v +
At F
4: r. .++1 r
Atv, k+
5: Compute pt x+1 using Eq. 13
6: Compute F. k+ (p X k using Eq. (1)
b+ At F *+4 0,
7: 0
k++
repeated for each time step, l E {0, Al,... Courant-Friedrichs-
(.....T), 7). where the time step. At, is chosen according to the (CFL.) condition, At 0.4h/c. This algorithm has the following physical interpretation: it prevents spatial information transfer through the code at a rate greater than the local speed of sound (small in the almost incompressible case considered this manuscript).
B APPENDIX: METHODS
B.0.1 FORWARD SENSITIVITY ANALYSIS
general, the loss functions this work can be defined as
L(X,0) = /
TT ( 0, d).
(22)
16

---

Under review as a conference paper at ICLR 2022
The forward SA (FSA) approach simultaneously integrates the state variables along with their sen- sitivities (with respect to parameters) forward in time to compute the gradient of L;
de dx // 0x
0, t)d..
(23)
Where, through using the chain rule, we see that the sensitivities the state variables with respect to the model parameters Xd are required to compute the gradient of the loss. Assuming that the initial conditions of the state variables do not depend on the parameters, then 0 (0) = O.
Now, define the sensitivities S dxx/d0% Then, from Eq. (8) we derive
dsp
dF. (X)),
(24)
d0a
then resulting, after applying the chain rule, in
dsa
OF (x(() S. OF (x 0)
(25)
di
0x:
80a
Since the initial condition X(0) does not depend on 0, then (0) = 0. Now, computing the gradient of the loss function reduces to solving forward time Initial Value Problem (IVP) by integrating simultancously the state variables X; defined in the main text, and sensitivities Sy, defined in Eq. (25). However, this means that the computational cost of the FSA method scales with (N)), where k is the number of parameters, which can be large when using NN's. With S. known, the parameters can then be learned (estimated) through an iterative first order gradient based approach (such as by updating using adaptive moment estimation (Adam) (Kingma & Ba 2017)) where the gradient is computed in Eq. (23)
In order integrate Eq. 25 the gradient of F; with respect to the parameters, both OF, (T),0))00)))))) 0)/00* and the Jacobian matrix, (0F;(X(7) 0) /0x, need to be computed. In this work, this done with mixed mode approach. OF (X 0) /0/8", with F. RR - R2d, computed with AD (the choice of forward reverse mode is determined by the dimension of the input and output space), where 1 is the number of parameters and d is the dimension. For example, >> 2d (as is the case when NNS are used), reverse mode AD is more efficient than forward mode (Ma et al. 2021). The Jacobian matrix is computed and obtained through mixing symbolic differ- entiation packages (or analytically deriving by hand), as well as mixing AD. For example, when there are NNs used for the parameterization of the right hand side, then according to expression for the Jacobian from the main text, AD derivatives will need to be computed on different functions each with potentially different dimensions of input and output space. With F, FX)) ) Rd + R2dd forward mode AD is more efficient. As seen the subsection 4.T, some of the models we train require neighborhood search list algorithms, which this case AD is only applied to derivative within the Jacobian, especially whenever summation over particles is done. Thus, a mixed mode approach is taken by combining both forward and reverse mode (depending prove the efficiency over just applying one or the other. The AD packages used in this work were both ForwardDiff.jl (Revels al., 2016) for forward mode and Zygote.jl (Innes 2018) for reverse mode.
dimensions) to im-
[0]2
I
[0]2 12 ( 0x,
OF, (x))((
OFF
OF
OFT
OF
(26)
0x;
d..
Ox7
Ox.
OF
OF
OF
OF
Or. 0x7
Ox
Ox
B.0.2 ADIOINT METHOD
This section provides an outline of the Adjoint SA (ASA) method used in this work (and can be found in (Bradley) 2019) (Donello et al. 2020)). Again, the goal is to compute the gradient of the loss function. This is continuous time dependent formulation, where the goal is minimize
17

---

Under review as a conference paper at ICLR 2022
loss function L(X (6, t), 0) which is integrated over time, ((X,0) S (XX 0,tdt subject to
the physical structure constraints (ODE or PDE), ( X, 0.t. = 0. and the dependence of the initial condition, g (0), 0) 0, parameters. Here, H is the explicit ODE in the standard form obtained through the SPH discretization of the PDE flow equations Eq. (8).
H X 0, t) x
(t),
(27)
gradient based optimization algorithm requires that the gradient of the loss function,
dg dgT(X,0)= (xx( 0 (x (0, 0g V X 0,(),
be computed. The main difference in the FSA and ASA approach is that in the ASA calculating de X is not required (which avoids integrating the additional & ODEs in FSA). Instead, the adjoint method develops a second ODE (size of which is independent of k) in the adjoint variable À as function of time (which is then integrated backwards in time).
The following provides a Lagrange multiplier approach to deriving this ODE in A. First define
L f= C= (V (=/_(v(X_0,0)) (V(X,0,))
0, (X, 1²² g ())
(28)
where À and are the Lagrange multipliers. Now, since H and zero everywhere, we may choose the values of and j arbitrarily. Then, VoL = Vol, resulting
Vol = 1."
( (Ox X + vo + X (0 Hdo 2x Hdo + 00H1) di
+ (8x10)dde 0øg).
Integrating the equation by parts, and elimitating, do we arrive at
Vol = 1/" [[xxV (0x d,0g H) x H) do + +0. + 2T dt
H +--XTa H gt) do do 009
",
Since the choice of, X7, and, μ. arbitrary, we set, X (T) = 0. and, -M (AT ())))))))))) (o))) in order avoid needing compute, de X (T), and thus canceling the second the last term in the latest (inline) expression. Now, assuming that the initial values of the state variables, X (0), do not depend on the parameters, we derive. dg (0) 0, and, g O. Then, we use loss function that does not depend on 0 explicitly, so that, 0g L = 0. And finally, we can avoid computing X at all other times > 0 by setting
xx + XT (0x 0,0x H)
H²
(29)
The resulting equation for the time derivative of À can be re-stated as the following adjoint ODE
OF
X = 0x XT
xt,))=0)
(30)
0X
where we also used that according to Eq. Combining all of the above, see that the simplified equation for the gradient of the loss function
27 0x H = - and 0 H = I2a-
is
Vol = VoL J
XToo Hdt.
Which, according to Eq.
27 becomes
1d1 /1 8F
VoL
X
00 dt.
(31)
18

---

Under review as a conference paper at ICLR 2022
Therefore, in order compute the gradient, the IVP expression (30) needs to be integrated back- wards in time for A(t), starting from A(ty) = Similar to the FSA formulation found above, both F and F,FF are computed with mixture of forward and reverse mode AD tools, depending the dimension of the input and output dimension of the functions to be differentiated.
Although the adjoint method has computational cost that is independent of the number of parame- ters, it requires more memory as the adjoint Eq. (30) must be solved backward in time (and requires the forward solution be stored). Notice that this workflow is not adequate for problems where real-time sensitivities are needed, as discussed in (Donello et al. 2020). This usually means that for (110)) FSA more efficient but when & > 100, the situation reverses and ASA becomes more efficient. See (Ma et al. (222]) for more details on the difference in efficiencies of mixing FSA (or ASA) with forward and reverse mode AD.
B. LOSS FUNCTIONS
In the previous subsections, we introduced the FSA and ASA methods for computing the gradient of the loss function, which used in our gradient based optimization learning algorithms. this section we construct three different loss functions; trajectory based, field based, and statistical based. The trajectory based loss function tries to minimize the difference in trajectories of the predicted particles from that of the ground truth trajectories. The field based loss function tries to minimize the difference in the velocity fields of the predicted and ground truth data (more fitting for the goal of SPH: approximate the field quantities). The statistical based loss function tries minimize the difference of underlying probability distributions obtained from the predicted and ground truth data. We explored all three loss functions and combinations of each this work. However, since our overall goal involves learning SPH models for turbulence applications, it is the underlying statistical description we want our models to learn and generalize to. This discussed further in section where we compare the hierarchy of models that is constructed in the next subsection 4.1
B.1.1 TRAJECTORY BASED LOSS
A simple loss function to consider is the Mean Squared Error (MSE) of the difference particles and velocity:
1 MSE = N X - 11
where X and X the ground truth and predicted particles states (both position and velocity) over the N particles respectively.
B.1.2 FIELD BASED LOSS
The field based loss function tries to minimize the difference of the large scale structures found in the fields of which the particles used to approximate. This is obtained by interpolating the velocities of the particles onto a mesh as seen in the following.
L,(0)
N,
I' -11,
where
N₂ my L Wo, (1lr// 1))
V
uses the same SPH smoothing approximation to interpolate the particle velocity onto a predefined mesh r- (with N, grid points), and V the velocity field prediction extracted from SPH samples (multi-particles temporal snapshots). The motivation for using field based loss function comes from the desire to reconstruct the underlying field variables (which is what SPH approximating for).
19

---

Under review as a conference paper at ICLR 2022
B.1.3 STATISTICAL BASED LOSS
In order to approach turbulent flows, one can well established statistical tools/objects. In this direction consider the integrated Kullback-Leibler divergence (KL) as a loss function (that will later be used in comparing hierarchy of models).
L( - f/
Ggt (gr()) (n)(t,,)) log prea,,l, ( (1)
dxdt,
(32)
where (t), and 2pr (t) represent single particle statistical object over all particles (such as the single particle dispersion statistics or velocity increment) over time of the ground truth and predicted data respectively. Here G(t, (f), .>) is the continuous probability distribution (in x:) constructed from data (t)) using Kernel Density Estimation (KDE). KDE a classical unsupervised statistical learning method in order to obtain smooth and differentiable distributions from data (Chen 2017). with
G(T, ((T,,,I)= Nh Σ K ((==x)//h)
(33)
where K is the smoothing kernel (chosen to be the normalized Gaussian in this work). Note the similarity with the SPH smoothing process, which makes it natural choice in this setting.
The integration in time becomes summation over discrete time steps once numerical integration over time is done (see|subsection A.4) in simulating the particles and performing SA (seelsection 4). To simplify the resulting expression, let us fix the time then integration over time can be carried last (through applying Leibniz integration rule). One derives
Gye(T, Z)) dx.
I(0,r)=
Gyr(T, Zyt, x)) log
Gpred (T, 2pr,
Combining all the preceding formulas we derive
ΣN N2 K
Zgt
N1
N1hu
h1
L(O, =
K
ggt
log
da.
N
N1h1
h1
Eppr
K
N2h2
Li=1
h2
For notational convenience, define the integrand as kl
-i=1 K
N1
99 -
Nih1
h1
kl(T, Fgt: Zpr,
K
2gt
x
log
N1
h11
1
N2
No
K
Z
Nzh2
2i=1
112
and
/
100
L(0)= =
ki(TT Zgt: Zpr: 2 )dz
so that
kl(TT Eyt; Epr, x)) IR TRRN
x
R N2
R
In order to perform gradient descent, this loss function will need be differentiated with respect to each parameter, and through the chain rule we will need the derivative of the state space variables with respect the parameters (found through FSA or ASA mixed with AD). Some of the follow- ing derivatives are done through automatic differentiation, however, in order to not require taking derivatives through the non-differentiable neighborhood list algorithm, some of the derivative need to be done analytically (or at least symbolically).
aL
OKl(T, Zgt, Zpr, x))
Økl(T, T, Egts Zpr x) Ozpr dx.
+
00a
Ozg
Ozpr
00a
20

---

Under review as a conference paper at ICLR 2022
where
Mi
K
ge x
Okl
1
gt
0
Nih1
1=
h1
K
log
i Eprs
N1h1
h1
Ozj
1
App
N2
K
N2h2
-==1
h12
Apr
K'
N1
Nah2
h2
K
gt
x
h2
N1h1
h1
im1
Li=1
K
App 1
N2l12
h/2
Zpr: 2
K'
1
N2ht2
h2
Gul
Gu(Zzt
hv
Gpred
(Zpry 2))
Using the velocity increment in both u and e this can be simplified:
aL a0a
N₂
Ökl, Ouh. Ozu
Okl, d da Ozy, 00%
00
where xp. = (T) ui (0). Simplifying,
0L
Our
000
Okly O2n
Ov 06
n Okly
Σ
dx
dx
00%
00
1=1
Out
Ov
Ökl, 02p
where
and
are computed using our mixed mode approach (in section 4
and
00
Ogo
Okl,
are computed with forward mode AD, and the integrals over x computed using numerical
Ozh integration.
21

---

Under review as a conference paper at ICLR 2022
B.2 LEARNING ALGORITHM
In what follows, we combine all of the computational tools and techniques introduced so far to outline the mixed mode learning algorithm.
Algorithm 1: Mixed Mode Learning Algorithm
Given Ground Truth: {X (lo), (4)) ... (())) SPH data Select model: d, (t) = (t), 0) from hierarchy
Select SA method (Forward Adjoint)
Choose Loss methods (Trajectory, Field, or Probabilistic) Choose optimizer: e.g Sprop ADAM
for & E n1, do
Prediction step: Verlet integration of model X Verlet(FF, 0x, ((t)), to,t/) ky)
OF
OF
Simultaneously computes
with mixed mode AD.
ax,
00
Forward SA then
Simultancously integrate system of ODEs for sensitivities S.;
0F (x, 0x) S. +
OF (X(((()). 0k)
d.S.
0x,
00,
else if Adjoint SA then
Integrate d
backwards in time (store forward solve);
d = 0x V X OF
X(t,) = 0
0X
end
Obtain predicted probability distribution & = KDE(2)
Compute VL (depends on SA method)
update using optimizer;
0 h++1 Opt(Ox)
end
Result: Estimates θ so that Lagrangian model is fitted SPH data with locally minimized loss
C ADDITIONAL RESULTS
22

---

Under review as a conference paper at ICLR 2022
PIMLSPII: 1096, timeestepe 10
PIML SPH Loss: N 4096
1000
10-3
10- 10
Pananeter
Loss
10
10-
2.5
25
10-
2000 4000
Iterations
6000 8000
2000
Iterations
4000
6000
8000
Figure 6: Solving inverse problem for the fully physics informed model 3D SPH flow with deterministic external forcing on 4096 particles over 4 SPH parameters. The solid lines show the SPH model parameters converging the dashed lines representing the ground truth parameters. Here the field based loss function (see ssussection 4.4) is used and is converging up to the order of machine precision (as Float64 and the MSE of field difference was used).
WC SPH steps
WC_SPH AV: 4096, time steps
12
10"
10
10
Loss
10
30
20
30
2000
3000
4000
0903
Iterations
iterations
Figure 7: Solving the inverse problems for 3D SPH flow with deterministic external forcing on 4096 particles over parameters (including the rate of energy injection). These plots shows the ability of the mixed mode method to learn the physical parameters of the Physics Informed model (Equation 7 The solid lines show the SPH parameters (initially chosen be uniformly distributed about (0.1)) converged to the dashed lines representing the ground truth parameters. Here the field based loss function (see subsection 4.4) is used with our mixed mode gradient based learning method (see section 4); we see the field based loss (b, d) converging up to the order of machine precision (as Float64 and the MSE of field was used).
23

---

Under review as a conference paper at ICLR 2022
Loaming
X Pa
p)
P
Figure 8: Using a Neural Network to approximate the equation of state (see subsection 4.1 using Lyz + Ly loss function. see that the Equation of State (EoS) P(p) is well approximated on the Pressure deviates once densities go
the domain of densities that is seen in training, however, far beyond range seen in training. Here Pot is the ground truth EoS and Pan is the neural network approximation of the EoS. We note that the underlying ground truth data has density variations within 1% of mean density, so the NN sufficiently approximates the EoS over a larger domain of densities seen in training.
24

---

Under review as a conference paper at ICLR 2022
Training set
Generalization over time
Truth
Physinf
EOSNN
GradP
Rotlav
NNSSm
NODE
Figure 9: Snapshots of flows comparing the ground truth the learned models, where learning is done with the LKI. + Lf loss function. NODE represents the least informed and Phys Inf is the most informed model. The main takeaway here is that even though the training occurring on shortest (physically relevant) time scale 1x (which corresponds to the time of neighboring particles to double their relative separation) the more physics informed models are able to generalize to much longer time scales all the way to the longest (physically relevant) time scale, beddy (which is the tumovvr time scale of the largest eddy of the flow).
25

---

Under review as a conference paper at ICLR 2022
Gonrslization Firee Over t
Rotutiona- - Inveriamm: Emrux
(PFRX)
Method
Method
(a) Generalization over time
(b) Testing Rotational Invariance
Figure 10: 2D generalizability where the models are trained and tested using the statistical based loss function. We also that the rotational invariance (and therefore conservation of angular momentum) is preserved as more physical struture is included in the model
26


Reviews:
[{"review_id": "fDzP31R5S9i", "reviewer": "Reviewer_hB3G", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "Strengths:\n\n- The identified problem lies within a very active area of research which is of huge importance in industry.\n\n- The paper is presented in a clear, logical manner with very good supporting supplementary material.\n\n- The methodology is valid: NNs are suitable for regression problems; automatic differentiation and adjoint-based sensitivity analysis are both used for gradient-based optimisation; the proposed loss functions are appropriate.\n\nWeaknesses:\n\n- The most-physics-informed model in the hierarchy is not NN-based, but has a few physically-interpretable parameters to be inferred, which can be done with standard gradient-based optimisation techniques. What are the benefits of using the less-informed, NN-based models over this one? Will the NN-based models perform better on real data? If not, it is difficult to argue their use case.\n\n- It would be worth discussing the relative strengths and weaknesses of the various models, for example by weighing the relative importance of model flexibility offered by using NNs versus model accuracy (or generalization/other metric).\n\n- The quality of the figures could be improved: the axes labels, tick labels and legends are too small to read in every figure, making it difficult to interpret results.\n\n- The sections describing forward, adjoint-based and mixed-mode sensitivity analysis do not add much insight to the paper (as they are not novel, according to the authors), nor are the efficiency benefits over standard automatic differentiation in these models evaluated. Would it not suffice to mention that forward and adjoint-based sensitivity analysis were used to increase the efficiency of gradient calculations (and include details in the appendices as necessary)? \n\n- The MSE loss is introduced but does not appear to have been used. Consider removing it.\n", "summary_of_the_paper": "The authors propose a hierarchy of physics-informed, mostly NN-based models of turbulent fluid flow. The NNs are trained to learn various unknown quantities (accelerations or pressures) of the turbulent system’s governing equations (which are derived in the smoothed particle hydrodynamics (SPH) framework) from a synthetic dataset, generated with a fully physics-informed model with known parameters. Automatic differentiation and adjoint-based sensitivity analysis are used to calculate the gradients of the loss functions efficiently depending on the relative sizes of the input and output spaces. The authors compare each model in the hierarchy in terms of accuracy, generalization, interpretability and data efficiency, and show that the more a model is physics-informed, the better it performs on these benchmarks.", "summary_of_the_review": "Although the paper is technically sound and the results are correct, the use of neural network-based models is not well motivated, given the fully physics-informed model (which outperforms all of the NN-based models) of the system does not require NNs (only standard gradient descent) to find its parameters. Comparing the relative strengths and weaknesses of the different models, and arguing the use-case for the NN-based models would greatly improve the paper. The paper demonstrates the methodology on synthetic data but would be more convincing with an example of it applied to real data.", "strength_and_weaknesses": "", "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": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "ZYE3s_5L5C5", "reviewer": "Reviewer_mpMw", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "Overall, the paper is well written, providing a very good introduction and overview of the modeling tasks. However, I have some concerns regarding the novelty of the work (wrt the methods), that are reflected in my score \n\nstrengths:\n- a good overview of the problem and motivation of the study\n- SA autodiff\n- reproducible\n\nweaknesses:\n- the contributions seem minimal \n- lack discussion and evaluation on some aspects of the considered approaches wrt. efficiency, scalability, robustness to noise etc.\n\nQuestions and comments:\n\n- Can you provide more the details about the fully-physics-informed model, and how are the parameters estimated? Are the same parameters used for each particle when generating the data? How would variances/noise in these parameters affect the performance of the models (I don't expect much change in the fully-physical model, but I'm curious about the other)\n\n- It seems that adding the artificial viscosity $\\Pi$ improves the generalizability of the models. How was it obtained, were the ground-truth values of $\\alpha, \\beta$, and $c$ used? \n\n- What is the compute time for each of the NN-based models and how well they scale on problems with more particles and longer trajectories? In this context, in the case of NODEs, how much is m, and how this affects the performance (larger vs smaller neighborhoods)? How much was l?\n\n- Are the same models reported in Fig4 used in the generalizability analysis? What is the performance of the models wrt (combinations of) different loss functions? It seems that the trajectory loss wasn't evaluated, is there a reason for this (other than potential over-fitting)? \n\n- As stated, the mode of the mixed-mode autodiff is chosen a priori for each model separately. But it seems that in all the experiments only FSA was used? How does ASA compare in such cases wrt efficiency?  \n\n- Can you elaborate more on the peaks in the loss in Fig.3?\n\nminor:\n- Can you elaborate more on the relation between the proposed PIML models and the ones outlined in the related work, that also relay on prior knowledge but attempt at learning parameters (eg. de Anda-Suarez et al 2018)\n\n- I would suggest moving the 'Why PIML challenges' paragraph from the Introduction to the discussion at the end (maybe?). While I enjoyed reading it, it's placement seems to break a bit the flow.  ", "summary_of_the_paper": "The paper explores and discusses the effects of incorporating prior domain-knowledge for modeling fluid dynamics with neural networks. The authors focus on smoothed particle hydrodynamics (SPH), a Lagrangian approach for modeling flows, and investigate different alternatives for incorporating physics-based inductive biases when constructing the models. To this end, the paper presents five different alternatives of such models that model particle velocities, with varying levels of prior knowledge: from fully data-driven Neural ODEs to fully physics-informed models with known structure (but unknown physical parameters). The authors also provide details on training such models using sensitivity analysis-based forward and adjoint methods, that are used interchangeably w.r.t the type of the model. Several experiments on synthetic data, in both transductive and inductive setting, show that model-variants that rely more on prior domain-knowledge are more accurate and capable of conserving quantities such as linear and angular momentum. ", "summary_of_the_review": "Overall, the paper is well written, providing a very good introduction and overview of the modeling tasks, which further motivate the importance of using prior knowledge in such scenarios. That said, while the study seems correctly executed, the overall technical contributions of the work seem marginally significant. In particular, in a more general sense, the benefits of knowing particular characteristics of the dynamics (or the whole structure) and using them for modeling dynamics are known. However, in many real-world scenarios this is not the case, which leads to different approaches for data-driven modeling. The authors set-out to evaluate this trade-off (in a synthetic setting), but do not provide deeper discussion between the model-variants across different aspects (such as efficiency, accuracy, memory, scalability, robustness (to noise) etc) which might be of more interest for the ICLR audience.", "strength_and_weaknesses": "", "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": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "I6lE47rrQ7K", "reviewer": "Reviewer_jmEb", "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": "Based on the smoothed particle hydrodynamics (SPH) method, this paper proposed and explored a fair amount of models with different levels of prior knowledge embedded.\nWhen evaluating the learned results, it's good that this paper also considers different kinds of physical properties and the model's generalization ability, other than just comparing the trajectory.\n\nHowever, in the experiment parts, there should be some comparisons with other existing learning-based fluid simulators, especially the methods that also use the Lagrangian view, to show how well the proposed physics-informed model performs compared with others. For better comparison, it's also helpful to show some results of the training and inference time used by different models.\n\n", "summary_of_the_paper": "This paper introduces a hierarchy of parameterized models based on one of the Lagrangian simulation methods, Smoothed particle hydrodynamics (SPH). \nSeveral models with different levels of physics information embedded are proposed. \nWhen training these models, sensitivity analysis and automatic differentiation are used to compute the gradients, and different loss functions (trajectory-based, field-based, and KL divergence) are used.\nIn the experiments, this paper evaluates the ability of the parameterized SPH models to learn the system's properties from data and to generalize to different environments.\n", "summary_of_the_review": "Although a certain amount of experiments is provided in this paper, it lacks some comparison results to show how well the proposed model performs compared with other existing methods. \nAdding SPH-related physics information to the model architectures can be a good direction to improve the data-driven simulator, but the models used in this paper seem to be simply replacing some parts in the classical SPH solver, which makes me have some concern about the novelty. ", "strength_and_weaknesses": "", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "FraZIy7cFkJ", "reviewer": "Reviewer_sFNH", "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": "The paper claims two contributions: (1) presenting a new physics-informed learning method, i.e. the models in 4.1, mixed SA+AD differentiation, combining MSE+Statistics-based loss, and (2) analyzing these methods, i.e. solving \"inverse problems\" and the finding that adding physics knowledge increases accuracy and generalization.\n\nOf the models in 4.1, the first one is Neural ODE with a per-node MLP, and the last is the actual PDE that is being learned, with 4 learnable system parameters. The models in-between successively add parts from the PDE into the model (summation, viscosity, kernel, etc.). The most interesting model version is the rotationally equivariant NN (though not exactly novel, given the large amount of papers on SE(3) equivariance in the last few years)\nInjecting physical knowledge into learned models is an interesting topic to research, but I was a bit puzzled by this particular approach.\nThe premise generally is that either we don't know the exact model equation, we have partial or noisy observations, or that ML components can give us other advantages over directly solving the PDE (i.e. faster inference, larger time-steps, differentiability, etc.). Otherwise there seems little point in attempting to use ML in the first place. Neither of these seem to be studied here-- the methods are trained and tested on the exact same PDE that is hard-coded as \"physical knowledge\", not even noise is added to the data. It is completely unsurprising that the methods which are closest to the PDE that generated the data perform and generalize best; and that the model which have the least learnable components train faster. Beyond this finding, there didn't seem to be any recommendations or pros/cons for the different methods in 4.1.\n\nIt's also worth pointing out that there are quite a few papers which do learn from SPH simulation, and already study a wide range of prior knowledge, from strong (the cited Ladicky et al.) to weaker (the cited Ummenhofer et al. and many GNN-based methods, such as Sanchez et al. \"Learning complex physics using Graphnets\"). It would be a good idea to study this field in context of these established methods as well, and ideally perform comparisons.\n\nThe \"inverse problem\" solved by the fully-physics-informed method is really just the equivalent of training the model-- i.e. fitting the 4 scalar models parameter to the data.\n\nI'm not sure I fully understood the implication of mixing SA and AD, and I'm not sure why this is necessary in the first place; the 'NODE' model is quite literally Neural ODE, so training the model in the same way as Neural ODE would be the obvious choice. Or, alternatively, directly supervising on dri/dt, dvi/dt. (the hard-coded model part should be trivial to implement in a differentiable manner).\n\nThe statistics-based loss is an interesting addition-- unfortunately the paper didn't study its effect, or explain why this is necessary when other methods can learn SPH simulations with a plain MSE loss.", "summary_of_the_paper": "The paper studies several methods for predicting SPH simulations of forced turbulence, in a Neural-ODE style setup. The methods range from general (MLP) to physics-based (i.e. the actual underlying PDE, with a few learnable parameters).\nThe main finding is that the more physics knowledge is provided, the more accurate and data efficient the learned model becomes.", "summary_of_the_review": "I don't think the paper in its current form contributes much new insight for the ML community. There isn't much technical innovation in the presented model, and the main finding is quite obvious-- under the conditions studied there really isn't much of a case for using ML methods at all.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "-z9LugCg6Mh", "reviewer": "Reviewer_Xkyx", "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": "I’d like to thank the authors for a well composed article. It studies an important direction of applying machine learning to scientific computation and does so in a systematic manner. The main decision is based on strengths and weaknesses expanded below.\n\nStrengths:\n* Authors consider a wide range of models;\n* Authors consider a relevant set of generalization properties (longer time, different Re)\n* Authors provide detailed account of the methods employed to carry out the research + source code\n\nWeaknesses:\n* The main conclusion of the paper which is that including additional physics structure improves performance is assessed in a relatively simple setting and rather expected\n* The generalization tests do not fully expose limitations of the learned models in the SPH approach, making it hard to assess what big challenges remain/solved using this formulation.\n* The use of additional loss functions are of exploratory nature.\n* Limited assessment of strengths and weaknesses of the learned Lagrangian vs learned Eulerian approaches.\nSystematic comparison of learned simulators is a very challenging task, as the phase space is extremely large. In the present work the models are mostly compared based on loss functions derived for training, as well as statistical distribution of the velocity increments. The losses compare expectedly, with models that are computationally closer to the reference (ground truth simulator) performing better. A better understanding of the structural differences between the errors made by different models could increase the significance of this result.\n\nGeneralization tests did not go into details of how models fail under the distribution shift. What is a particular challenge that less structured models struggle with? At which point/type of distribution shift even more structured models start to fall flat?\n\nThe additional (to my understanding novel) loss function, such as field reconstruction loss are interesting, but their use was not scrutinized. \n\nFinally the last important point brought up by the paper is the use of learned SPH methods, which is comparatively underrepresented in the literature of learned fluid dynamics models. To generate more interest in this approach a performance comparison between two methods or a clear scoping of what are the central challenges well addressed by the SPH method would be valuable.\n", "summary_of_the_paper": "The paper develops a learned fluid dynamics simulator based on smoothed particle hydrodynamics approach. Authors consider a range of models that parameterize the underlying equations of motion for individual degrees of freedom at different levels of granularity (originating from the variation in the amount of physics priors included in the model). A mixed mode AD is used to perform gradient-based optimization of model parameters to minimize discrepancy between reference ground truth data. Authors use a combination of three loss components (particle trajectory, field reconstruction consistency and statistical similarity).\n\nAuthors train and evaluate models on a compressible 3D turbulent flow with external forcing. A range of simulation parameters is considered during evaluation, assessing the ability of all models to produce accurate predictions over longer time horizons, as well as simulating more complex flows (different Reynolds numbers). The generalization performance is assessed based on the sum of field reconstruction error and statistical consistency loss. Additional comparison is provided on symmetry properties of learned models, as well as visual comparison of model predictions when trained on very short time sequence;\n", "summary_of_the_review": "While the paper reads well and presents results carefully, I found it hard to formulate a clear take away from the paper. This aspect is reflected in the main review section and motivates my rating of \"not good enough\".", "strength_and_weaknesses": "", "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": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper explores and discusses the effects of incorporating prior domain-knowledge for modeling fluid dynamics with neural networks, with a focus on smoothed particle hydrodynamics. Reviewers agree that the contributions are modest, and that they are not well presented with respect to issues of efficiency, scalability, robustness, etc.  More work need to be done to make it useful to the community.  Many of directions to improve the paper were in reviewer comments.
Decision: Reject