Enforcing physics-based algebraic constraints for inference of PDE models on unstructured grids
Authors: Valerii Iakovlev
Markus Heinonen
Harri Lähdesmäki
Venue: ICLR 2022 Submitted (2022)

Abstract: Data-driven neural network models have recently shown great success in modelling and learning complex PDE systems. Several works have proposed approaches to include specific physics-based constraints to avoid unrealistic modelling outcomes. While previous works focused on specific constraints and uniform spatial grids, we propose a novel approach for enforcing general pointwise, differential and integral constraints on unstructured spatial grids. The method is based on representing a black-box PDE model's output in terms of a function approximation and enforcing constraints directly on that function. We demonstrate applicability of our approach in learning PDE-driven systems and generating spatial fields with GANs, both on free-form spatial and temporal domains, and show how both kinds of models benefit from incorporation of physics-based constraints.

Under review as a conference paper at ICLR 2022
ENFORCING PHYSICS-BASED ALGEBRAIC CON- STRAINTS FOR INFERENCE OF PDE MODELS ON UNSTRUCTURED GRIDS
Anonymous authors Paper under double-blind review
ABSTRACT
The lack of relevant physical constraints in data-driven models of physical systems, such as neural network parameterized partial differential equations (PDEs), might lead to unrealistic modeling outcomes. A majority of approaches to solving this problem are based on forcing model satisfy a set of equations representing physical constraints. Currently available approaches can enforce very limited set of constraints and are applicable only uniform spatial grids. We propose a method for enforcing general pointwise, differential and integral constraints on unstructured spatial grids. Our method based on representing a model's 's output terms of a function approximation and enforcing constraints on that approximation. We demonstrate wide applicability and strong performance of our approach data data-driven learning of dynamical PDE systems and distributions of physical fields.
INTRODUCTION
Multiple works have shown the capability of neural networks solve complex physical problems and learn the behavior of physical systems from data. Examples include learning and solving ordinary differential equations (ODEs) [6], partial differential equations (PDEs) [28; 20] and rigid body dynamics [31; 5]. Purely data-driven models are typically not forced to satisfy physical constraints of the system that generated the data. This might lead to unrealistic predictions that violate some known properties of the underlying physical system.
Incorporation of relevant constraints allows to make a better use of the available data and makes predictions more physically plausible. The field dealing with poostre-conttrained learning is diverse and offers many approaches to adding constraints to models. We refer the reader to many reviews for details [30; 3: 36; 19]. The approach consider in this work is based forcing model to satisfy algebraic constraints represented by set of equalities and inequalities. This is the most commonly used upproach which allows represent wide range of constraints and has been shown work well in many cases [18; 17; 25]. However, while many constraints can be represented algebraically, it not always clear how evaluate and enforce them.
Currently available approaches to enforcing algebraic constraints limited to uniform grids and have a very narrow range of constraints they can enforce (e.g. only pointwise, or specific differential constraints), see Section 5 for details of related work. Such approaches to models based on convolutional neural networks (CNNs) but cannot be extended to recently developed models based graph neural networks (GNNs) [33; 27; 15] and other models working unstructured grids.
be readily applied
We propose a much more general method which allows enforce pointwise, differential and integral constraints on unstructured spatial grids and demonstrate its applicability in learning of PDE-driven dynamical systems and distributions of physical fields. The method is based on using a models's output at the nodes of grid to construct interpolant and applying constraints directly to that interpolant (Section 3).
Code and data will be made publicly available.

---

Under review as a conference paper at ICLR 2022
2 BACKGROUND
PDE-driven dynamical systems. Many physical systems can be described in terms of PDEs. Such systems defined bounded domain which they evolve over time. We consider continuous dynamical systems with state u(f, x) E Rº that evolves over time E I220 and spatial locations æ E Q C RD. For physical systems, D is typically limited to {1, 2,3 although our method will work with any value of D. We assume the system is governed by unknown PDE
@u(1, x) F(x, u(t, x), Va att Φ) ())
0t
which describes the temporal evolution of the system in terms of the locations a; state and its first and higher-order partial derivatives w.r.t. æ. The goal of data-driven PDE model is to learn the dynamics F' from data.
Data for learning F is collected by measuring the state of the system at observation locations (......., .N) over increasing time points (to. ........t.. ...)..).. This results in dataset (ft))...... y.tt)) where (t,)) = (u(ty, æ1).... u(t,, æx)) is a collection of observations. The dataset is used to train the model to predict (y(4), ( .......)))))))) starting from the initial state y(o)). Training is typically done by minimizing average loss between the model's predictions u(t) and the data y()).
PDE models differ restrictions they impose on time points (temporal grid) and observation locations (spatial grid). Some models require both grids to be uniform [23], other models relax these requirements and allow arbitrary spatial [27] and spatio-temporal grids [15]. We build our
algebraic constraints method using the model from [15] as the most general one. The model is based on application of the method of lines [32] to Equation which results into a system of ODEs
du(t.ar)) dt
F0(21, U1, UN(())
út) =
du(t..E
FolaN, UN, UN(())
which approximates the solution Equation at the observation locations x; using their neighboring points N(i), where (x)(() and (N(()) are the neighbors' positions and states respectively, and и; is u(t,æ;). The approximate solution converges to the true solution as N increases. The true dynamics F' is approximated by a parametric model F' whose parameters 0 learned by minimizing the difference between the model's 's predictions
u = u(0) + [' ('((rddr
(3)
and the data (t). The integral in Equation 3 is solved using a numerical ODE solver. In [15], the function Fg was represented by a graph neural network (GNN) which takes states and locations at an observation point and its neighboring points (ii).. The observation points are connected into grid using Delaunay triangulation which allows to naturally define N as set of points connected to the point i. However, F' can be represented by other models and a different neighbor selection criterion can be used. The model parameters 0 are learned by minimizing the MSE between (t) and u(t)
M
Laata
Σ
|[(())
(
(4)
M
The gradient of Laam w.r.t. 0 is evaluated using the adjoint method shown in [7].
Generative Adversarial Networks One of the tasks that we consider is learning distributions of physical fields. For that purpose we utilize generative adversarial networks (GANs). A GAN is generative model consisting of a generator and discriminator [12]. The generator, G, learns to transform a random variable ~ pz over a latent space Z to the data space y in such way that the discriminator, D, cannot tell the difference between samples generated by G and samples from the data distribution Pana. Both, and are learned by solving the following minimax problem
min max VGG, D) Ey [log D(Y)) + Ez~pz [log (1 D(G(Z)))
(5)
G
Solution of this problem exists and is unique with the optimal generator perfectly mimicking the data distribution [12].
2

---

Under review as a conference paper at ICLR 2022
3
METHODS
In this section we presents an approach to evaluating pointwise, differential and integral constraints on unstructured grids. Then, we demonstrate how this approach can be used to enforce arbitrary soft and linear hard constraints.
3.1 EVALUATING CONSTRAINTS ON UNSTRUCTURED GRIDS
We assume the data y() is available at observation points (#1..... x and that model makes predictions u(t) at these points. We assume the predictions to be evaluations of an unknown underlying function. Since the underlying function is unknown, we cannot impose constraints it directly. Instead, approximate it by interpolant (t, æ) and impose constraints on (f, æ) (Figure 1). The approximation is constructed from u(t) by placing a basis function at each æ; and representing us (1, as
and time points
(M)
uf(t,x) Σ a (1))
(6)
where dy is a scalar basis function at æj and aj E R'. The coeffi- cients (t) are obtained from u(t) (see Section 3.4).
U1
Next, we show how to evaluate constraints on uy æ)) using basic building blocks. To avoid cluttered notation, we consider equality constraints and assume u(t, x), ) x E R. Generalization to inequality constraints, vector fields and higher spatial dimensions is straight- forward.
U2
U4
X1
X2
X3
Xa
Figure I: Example of approx- imating unknown under- lying function (green) by an interpolant (red) constructed from observations i u1....,
Pointwise constraints. Consider points (Z1, (....... 2x) ) in Ω on which a pointwise constraint (t, 2;)) = should be evaluated. + representable terms of finite
Assume the function number of functions 7' (uy (t, 2.)) ) R + IR indexed by m. For example, should the constraint be h(uy) = 3u, + ut = 0, then we would define " (uy) = uy. (2 = u and h(u)) 3 (1f)) + r2(u)) 0. Then, h can be evaluated by evaluating each Tm as
(m,(u(t,z))) 7/m
(a,(t)o Ymm (i,, a(t)),
where a+t) (a) ax))) K-by-N matrix with elements Dú.j = (z(().). and 2.,, is the i'th row of Φ.
Differential constraints. Consider the same setup as before but now huuu (t, 2,)) 0 consists of differential operators and is representable in terms of a finite number of functions (m.
8xy
+ IR indexed by m, where the derivative order could be different for each m. For example, should the constraint be h(u,) 3uy Uf
du? 0. then we would define ?1 (uy) = uy. 7/ (u,) u² Then, J can be evaluated by evaluating each ym (1y (t.2))) using the generalization of the chain rule (Appendix A) which contains only types O2
dhe (1)(
and M(u)) 3 3 mi (ur) + (ur)
of terms. The first type of terms
d'm duy
d'am can be evaluated using automatic differentiation
..... du
while the second type of terms
can be evaluated as du, geu de: dz
Fu
N
(z,)
Σ
(t)
(000) a()).
(8)
Oz!
Oz!
where (0)
(z.) Oz?
Mixed partial derivatives can be handled in similar way (Appendix A).
3

---

Under review as a conference paper at ICLR 2022
Integral constraints. Consider the same setup before but with (t,
Ia (tu = 0, where the function r R - is representable in terms of functions I' Eny(uf (t, 4:)) ) IR + similarly to the pointwise constraints. Then, Jg Tugg(t, aa can be evalu- ated using numerical integration technique, e.g. midpoint rule, Gaussian quadrature or Monte-Carlo integration, as
J. ((uy(t, (t. d) = Σ ((,
where K is the number of integration points, μ4 are integration coefficients which depend on the grid and integration method, and (t, 2;)) is evaluated as in Equation 7.
3.2 SOFT CONSTRAINTS
Soft constraints are implemented by minimizing the following loss Tadata + Xr(h(ur)), where À E IR M and Laua is defined as in Equation 4. We r(h(uu)) KM Σ D== (tu z,)) for point- wise and differential constraints and r(h(ur)) मे Σ huuu x))² for integral constraints.
3.3 HARD CONSTRAINTS
Our method allows implement hard constraints by projecting the interpolant (t, x)) to a subset of functions which satisfy the required constraints. Namely, if uy(t, æ) does not satisfy constraints and h, is projected a subset of functions which satisfy the constraint by solving the following optimization problem
min 6 EV
|uu us 11
s.t.
Mū,)) = 0,
(10)
g(û S O,
where the projection is denoted by (t, x) and V& is spanned by the basis functions.
Using the basis representation uy (t, x) = E a, (t))):( and úy(t, x) = E (,(t)o((z) we
can rewrite the optimization problem (10) as
min B(E)ERN
(ox(1) - ((t)))((((()) ((())
(11)
s.t.
h(ār) 0,
g(A,)
where ((t) =
and ... Ja (1())))(())
To train the model end-to-end the problem (11) should be differentiable. Agrawal et. al. [1] proposed differentiable convex optimization which could be used in this case if the problem (11) could be expressed in a DPP-compliant way (see [1]). To that, we restrict ourselves to constraints that can be expressed an equality or inequality between AB((1) and b, where A is a constant matrix and is a constant vector. This formulation admits pointwise, differential and integral constraints on untransformed u.y The objective function is convex since its Hessian is positive-semideffitt i.e. for any R.
T S
Σ b,==1
(urs, 0ンو5)
(3bb)) L2 0.
(12)
This allows to solve the problem (11) and differentiate its solution (1) w.r.t. xx(1). The model parameters are found by minimizing the following loss function L + Xl Xeproo where À E IR and Laaa is defined as in Equation 4 but with u(t.) replaced by û(t,) = (ūy (t,, æ1)...,, set Lpooj = N. Σn" Σソ [Iuy((y, ax)) - ú æz)11 The second term makes the optimization procedure prefer models that predict uy close to the feasible set of the problem (11).
(N), We
We note that the proposed approach currently limited to small-scale problems due to existing computational bottlenecks in the implementation of differentiable convex optimization [1].
4

---

Under review as a conference paper at ICLR 2022
3.4 BASIS FUNCTIONS
Selecting appropriate basis is crucial for efficiency and applica- bility of the proposed method. Ideally, the basis should allow efficient construction of (y(4,a) from u(!), contain no tunable parameters, and lead to sparse matrices Φ, 414) and o We con- sider bases from two families: Lagrange basis functions and radial basis functions (RBFs).
Lagrange basis functions do not have tunable parameters and have compact support which leads sparse + (p()) and áb For the piecewise linear basis the interpolant (t, x) can be constructed directly from the predictions by setting ax(t) u(t). However, constructing (4, æt) for a higher order basis, e.g. piecewise servation points
Figure 1D and 2D piecewise linear basis functions (colored) and function built from them (grey). Black dots represent ob-
quadratic, requires the model make predictions not only the observation points, but also at some extra points where the data is not available. In Section 4 we demonstrate one approach to solving this problem. After extending the state u(t) by predictions at the extra nodes, the coefficients at can be evaluated similarly to the piecewise linear basis. In this work we use piecewise linear (PWL) and piecewise quadratic (PWQ) bases. Examples of PWL basis functions are shown in Figure 2.
Radial basis functions have wider range of properties. Some RBFs have tumable parameters, some don't. The matrices क Φ(9) and & evaluated with RBFs typically dense, but RBFs with compact support exist (e.g. bump function). The interpolant uy (t, æ) can be constructed by evaluating ax()) = K utt), where K- is the inverse of the interpolation matrix the given RBF and Ky .y = o æ zj [I), where @ is an RBF and T; zj observation locations. In this work use the cubic RBF basis i.e. <(r) = 33.
We use PyTorch [26] to handle sparse matrices and to evaluate K --tu(t) in a differentiable way.
4 EXPERIMENTS
In the following experiments we use the relative error between the data y() and model predictions u(t) defined as [ and consider only soft constraints. We present an experiment with hard constraints implemented as shown in Section 3.3 in Appendix D. Data generation is described in Appendix B. Training, testing and modeling details are in Appendix C. All experiments were run on single NVIDIA Quadro P5000 GPU. All errors bars represent standard deviation of the results over five random seeds.
[
4.1 REPLACING EXISTING METHODS
In this experiment we take existing models which incorporate physics-based constraints in training and replace their constraint enforcing approaches with ours. We consider two works. First, [37] which trains a GAN to produce divergence-free vector fields using zero-divergence constraint. Second, [10] which predicts warping fields driving the evolution sea surface temperature by observing snapshots of the temperature fields (see Appendix C for more details). Both models work on uniform grids which allows them to evaluate constraints using finite differences. For comparison, we replace finite differences with our method and observe how it changes the models' performance. In both cases we use the PWL basis.
time while enforcing gradient and divergence constraints on the warping
For [37] we track the mean divergence and discriminator loss. Results of the original approach are as follows: mean divergence 0.079 and discriminator loss 0.091. With our method the mean divergence was 0.014 and the discriminator loss was 0.088, Both approaches results in similar discriminator losses but approach produces smaller mean divergence (smaller is better). Our method increased the runtime per epoch by 6%.
For [10] we track the total, divergence and smoothness losses which, with the original approach, were 0.139, 8.4 10-5 and 1.51 10 4 respectively. With our approach the losses were 0. 139, 8.3 10-5 and 1.51 10-4, respectively. Both methods produce very similar results. Our method increased the runtime per epoch by 30%.
Overall, replacing existing constraint enforcing approaches by ours on data from uniform grids resulted in comparable model performance, except for runtime which was slightly increased.
5

---

Under review as a conference paper at ICLR 2022
Node frac. 0.5
Node frac. 0.75
Node frac. 1.0
Rel. erro
0.05
bed cubie
0.05
0.05
0.04
0.04
0.04
0.03
0.03
0.03
30
60
120
30
60
120
30
60
120
0.004
0.002
0.002
Constr. viol.
0.002
0.001
0,001
0.000
0.000
0.000
30
60
120
30
80
120
30
60
120
Train size
Train size
Train size
Figure 4: Effects of amount of data and grid sizes on relative errors and constraint violations for the Cahn-Hilliard equation. All results are for the test set. Constraint violations are evaluated as the mean absolute violation of the constraint, Jo ut x x) - C| over all simulations and time points. In most simulations = 0.5.
4.2 CAHN-HILLIARD EQUATION WITH AN INTEGRAL CONSTRAINT We start with the 1D) Cahn- equation
du Øt = (u(1=u))
- -w) -22 v
(13)
which is known to conserve the state u i.e. h(u) = In u(t, 2)dx C at all time points, where C - So u(0, x dz constant. This is an example of conservation law which are abundant in nature and are important class of constraints that data-driven models of physical systems should satisfy. Conservation laws can be expressed in differential and integral forms and this experiment demonstrates how the integral form can be enforced, The constraint is evaluated using the midpoint rule as shown the previous section with single "/1 being the identity function. We use PWI., PWQ and cubic RBF bases and compare the results to unconstrained model.
For training we use 30, 60 and 120 simulations while the test set consist of 60 simulations. Simulations in the training/test data last for 0.0015/0.000 5/0.0030 seconds and contain 50/100 uniformly spaced time points. The full spatial grid consists of 101 uniformly spaced nodes. We randomly sample 50%, 75% and 100% of the nodes and train/test on the resulting (irregular) spatial grid. Training and testing is done with identical spatial grids. An example of a spatial grid
0.0
0.2
0.4
0.6
0.8
1.0
Figure 3: ID spatial grid for the Cahn-Hilliard equation. with 50% of nodes is shown in Figure 3. We evaluate the constraint on a uniform grid with 200 nodes placed on top of the original grid.
To learn the dynamics of the system we use the model from [15] (Section 2). We found that using GNN produced poor results. For that reason perceptron (MLP) which updates the state of each node based on the states of all other nodes in the grid (results for GNN are in Appendix E). The MLP contains two hidden layers with Leaky ReLU nonlinearities. The number of hidden neurons was set to the number of nodes in the grid.
represented the function Fg with multiplayer
The coefficients (x()) for the PWL and cubic bases can be evaluated directly from the model predictions at the grid nodes. But the PWQ basis requires extra predictions to be available between the nodes. This is problematic since there is no data at these points to guide the model's predictions. To solve this problem we introduce a small MLP which is applied to consecutive pairs of nodes. The MLP takes the states at both nodes and the distance between them as the input and estimates the state at the midpoint between the two nodes. The MLP is trained jointly with the main model and uses only the constraint-related loss term during training.
For testing, we construct the interpolant (t, x)) using the thin plate spline basis (ó(r) r log r) and evaluate the constraint on that interpolant. This allows to make a fair comparison between the unconstrained model and different bases and avoid biasing or ovefitting bases used for training.

---

Under review as a conference paper at ICLR 2022
Figure 4 shows results of the experiment. We observe that changing the node fraction does not significantly affect the relative errors but has noticeable effect on constraint violations, especially for the unconstrained model. Constrained models tend to show similar or better performance than the unconstrained model. Among all bases, the cubic basis consistently results in lower relative errors and constraint violations. However, the simpler PWI. basis often performs on par with the cubic basis, especially on denser spatial grids. We also observe that coarsening of the grid increases the constraint violation gap between constrained and unconstrained models and that this gap seems to not close as we increase the amount of training data.
The PWQ basis performs rather poorly on fine grids which likely due to a suboptimal approach to evaluating the state at the extra nodes. A better approach could consider not only pairs of points but also larger neighborhoods. Nonetheless, the PWQ basis achieves good performance on coarse grids which shows that piecewise bases of higher order could potentially be used to enforce constraints. This will allow to scale grids with large number of nodes due to sparsity of the constraint matrices and efficient evaluation of a.
4.3 HRAT EQUATION WITH MONOTONICITY CONSTRAINT
We impose constraints on 2D system governed by the heat equation 0u U for which the generated initial conditions (ICs) are monotone in one direction. Since the ICs are monotone, the state u remains monotone at all time points as well. We enforce the monotonicity constraint as de
ăt
du 2 0. The
constraint evaluated as shown in the previous section with 71 being the identity function.
For training we use 15, 30 and 90 simulations while the test set consist of 120 simulations. Simulations in the training/test data last for 0.1/0.2 seconds and contain 21/41 uniformly spaced time points. The full spatial grid consists of 1087 nodes. We randomly sample 33%, 66% and 100% of the nodes and train/test on the resulting (irregular) spatial grid. Training and testing done with identical spatial grids. Spatial grid with 100% of nodes is shown in Figure 5. The constraint is evaluated at the nodes of a uniform 51 x 51 grid placed on top of the original grid.
Figure 5: Grid for the heat equation.
To learn the dynamics of the system we use the model from [15] directly with the messaging and aggrega- tion networks being MLPs with a single hidden layer consisting of 60 neurons with Tanh nonlinearities and
0.04 sec
0.1 sec
the input/output sizes of 4/40 and 41/1 respectively.
Data
During testing, we use predictions of the models to construct an interpolant (I, æ) using the thin plate spline basis and evaluate the constraint on that in- terpolant. This allows to make a fair comparison between the unconstrained model and different bases.
Unc.
Figure 7 shows results the experiment. We observe that changing the node fraction equally increases rel- ative errors of all models and has noticeable effect on constraint violations, especially for the unconstrained model. Constrained models tend to show slightly higher comparable relative errors but noticeably lower constraint violations than the unconstrained model. The cubic and PWL bases perform equally well in this case. Similarly to the experiment in the previous section, we observe that coarsening of the grid introduces larger constraint violation gap between constrained and unconstrained models and that this gap seems to not close as we increase the amount of training data.
Soft
Figure 6: Comparison of data with predic- tions of unconstrained and constrained mod- els trained on 30 simulations, full spatial grid and using PWL basis for the constrained model. The predictions are for test case.
Figure 6 shows qualitative difference between predictions of constrained and unconstrained models. can be noted that predictions from the constrained model have noticeably smoother contours thus making the field more monotonous in the horizontal direction.
7

---

Under review as a conference paper at ICLR 2022
Node frac. 0.33
Node frac. 0.66
Node frac. 1.0
0.03
0.03
0.03
Rei. error
0.02
0.02
0.02
0.01
0.01
0.01
15
30
90
15
30
90
15
30
90
0.0010
0.0010
0.0010
Constr. viol.
unc
I
cubic
pul
0.0005
0.0005
0.0005
15
30
90
15
30
90
15
30
90
Train size
Train size
Train size
Figure 7: Effects of amount of data and grid sizes on relative errors and constraint violations for the heat equation. Results are for the test set. Constraint violations are evaluated the mean absolute violation of the constraint duy 2 over all simulations and time points.
dx
4.4 LEARNING DISTRIBUTIONS OF PHYSICAL FIELDS
We demonstrate the effect of adding constraints GAN when learning distributions of physi- cal fields on unstructured grids. We use Wasser- stein GAN (WGAN) [2] as a more stable variant of a GAN. We use MLPs as a generator and discriminator. Unconstrained and constrained models are trained for 1.2M iterations. Con- straints are enabled only after 600k iterations. Constrained models are trained similarly to the unconstrained ones but with a modified genera- tor loss defined as LG A In (1 + Lc), where LG is the standard generator loss and Le is the constraint- based loss. We define T.c the mean value of hu))), where is a constraint evaluated the centroid of each cell in the grid.
02
Ia unc.
0.
0.6
1.2
-2
Iteration
Divergence
0.1
Lg tine
0.0 noration
500
-25
25
50
4.4.1 ZERO-DIVER DIVE-DIVERGENCE FIELDS
Laplacian
Divergence-Fie vector fields are often encoun- tered in solutions of fluid dynamics problems. The divergence-fir free constraint on vector field u(x, y) = (ur (I. y), u2(2. (x)) y))" is defined as h(u) = d0
Figure 8: First row: Lp. Le and histograms of divergences of samples from data, constrained and unconstrained WGANs. Second row: Lp, Le and histograms of Laplacians of samples from data, constrained and unconstrained WGANs. Con- straints are evaluated cell centroids.
au + Pug = O. The constraint en-
by
forced using the PWL basis. We generated dataset with 10k divergence-free fields on an unstructured grid with 1050 nodes (Figure 14)
and used a WGAN to learn distribution over such fields. Note that the generated fields are not entirely divergence-free but have small residual divergence due to discretization errors.
Figure 9a shows that there is clear difference in the quality of the samples generated by the unconstrained and constrained models. Samples from the constrained model are smoother and more similar to the data. Quantitative comparison of the samples presented in Figure 8 shows that the constrained model generates fields that have much lower constraint violation and divergence distribution very similar to that of the data.
4.4.2 ZERO-LAPLACIAN FIELDS
Fields with zero Laplacian represent solutions to some PDEs, for example the steady-state heat equation. The zero-Laplacian constraint scalar field u(z, y) is defined as h(u) The constraint is enforced using the cubic basis the PWL basis has zero second derivatives
0. ag ayy
8
a²
8

---

Under review as a conference paper at ICLR 2022
co
40
40
Unconstraieed
Constraina
Constranee
bb)
Figure 9: Magnitudes of random samples from the dataset, unconstrained and constrained WGANs. a) -ire-diveeeence fields, b) Capl-Lapl fields.
everywhere. We generated dataset with 10k Laplacian-free fields on an unstructured grid with 1050 nodes (Figure 14) and used WGAN to learn distribution over such fields. Note that the generated fields are not entirely Laplacian-free due to discretization errors.
Results of the experiment shown in Figures 9b and 8. Similarly the divergence free case, visual quality of the fields generated by the constrained model is significantly better than for the unconstrained model. Quantitative comparison of the samples presented in Figure shows that the constrained model generates fields that have much lower constraint violation and Laplacian distribution very similar that of the data.
RELATED WORK
Soft constraints. Soft constraints are widely used due to being relatively easy implement. Examples include lake temperature prediction [18; 16], traffic simulation [22], fluid and climate modeling [11; ; 10; 4], where constraints are evaluated pointwise or using finite differences.
Hard constraints. Approaches implementing hard constraints diverse and can be categorized as processing the output of an unconstrained model [4; 17; 24; 34] and designing a model that produces feasible predictions by default [23; 25; 14; 13; 8: 38].
Constrained PDE models Current approaches to enforcing soft [10; 11] and hard [21; 25; 17] constraints are limited to specific types of constraints and spatial grids. For example, [25; 17] implement only hard differential constraints and both are limited to uniform grids. Uniform grids allow to evaluate constraints efficiently e.g. using finite differences [10; 21; 25] or fast Fourier transform [17] but assuming that the data lies uniform grid might be limiting.
Constrained GANs Works such as [37; 17] showed how physics-basee based constraints benefit training and quality of the generated samples are also limited to uniform grids.
CONCLUSION
We presented a general approach to enforcing algebraic constraints on unstructured grids and showed how it can be used enforee soft and hard constraints. We demonstrated applicability of the approach to learning of PDE-driven dynamical systems and distributions of physical fields. We considered two families of basis functions for constructing the interpolant and showed how Lagrange basis functions of order higher than one can be used. Our method allows to drop the unrealistic assumption about uniformity of spatial grids and shows promising results various tasks.

---

Under review as a conference paper at ICLR 2022
REPRODUCIBILITY STATEMENT
All details required to reproduce the experiments are provided in Section 4 and Appendices. Code and data used to run the experiments will be made publicly available after the review process.
REFERENCES
[]] A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and Z. Kolter. Differentiable convex optimization layers. In Advances in Neural Information Processing Systems, 2019.
[2] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In Doina Precup and Yee Whye Teh (eds.), Proceedings of the 34th International Conference Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp- 223. PMLR, 06-11 Aug 2017. URL http:/ http::/tttttceeeeeeeeeoo /ttp///oooddnngg ml r. ww.pprsssvvv arjovsky17a htmm.
[3] Nathan Baker, Frank Alexander, Timo Bremer, Aric Hagberg, Yannis Kevrekidis, Habib Najm, Manish Parashar, Abani Patra, James Sethian, Stefan Wild, Karen Willcox, and Steven Lee. Workshop report basic research needs for scientific machine learning: Core technologies for artificial intelligence. 2 2019. doi: 10.2172/1478744 URL https: /ttp:://wwwwwwe isstii..... gov/
biblio/ /////1477744
[4] Tom Beucler, Michael Pritchard, Stephan Rasp, Jordan Ott, Pierre Baldi, and Pierre Gentine. Enforcing analytic constraints in neural networks emulating physical systems. Physical Review Letters, 126(9), Mar 2021. ISSN 1079-7114. doi: 10. 1033physrevleet.. wttp://eee....0988.... 126.098302. URL
http:// / dx. doi. org/ 10. 1103 / htw...1///yyseee/////// Phys 126. 098302.
[5] Michael Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum. A compositional object-based approach to learning physical dynamics. arXiv preprint arXiv: 1612.0031, 2016.
[6] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL ht http:///ooceeeeee/e//e//c////cat neurips cc pccppappe/2018/ pdf.
file/ 6.3e66938655 6bb dddedddddd 698
[7] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL ht /ttptppee/pppccoeeeeeeee ngs neur cc pccppppp////01/ pdf.
file/ (.338fffbllfffeeec (tt.ee6662222444699562
[8] Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks, 2020.
[9] Arka Daw, R. Quinn Thomas, Cayelan C. Carey, Jordan S. Read, Alison P. Appling, and Anuj Karpatne. Physics-guided architecture (pga) of neural networks for quantifying uncertainty in lake temperature modeling. Proceedings of the 2020 SIAM International Conference Data Mining, pp. 532-540, Jan 2020. doi: 10. 00.113771.77816111966666
URL http:
/ //dx. doi oww./////00 111./1..71111111111 1711117118778881 973366 60.
[10] Emmanuel de Bézenac, Arthur Pajot, and Patrick Gallinari. Deep learning for physical pro- cesses: incorporating prior scientific knowledge. Journal of Statistical Mechanics: Theory and Experiment, (019(12):124009, Dec 2019. ISSN 1742-5468. doi: 10.1188/1744-5568//b333333 URL http: dxx doi. oro/ 10. 1088 1742 2 5468 (0/56///33/ ab3195
[11] N. Benjamin Erichson, Michael Muehlebach, and Michael W. Mahoney. Physics-informed autoencoders for lyapunov-stable fluid flow prediction, 2019.
[12] Ian Goodfellow, Jean Pbadet-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sher- jil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q. Weinberger (eds.), Advances in Neural Information Processing Systems, volume 27. Curran Associates,
10

---

Under review as a conference paper at ICLR 2022
Inc., 2014. URL
[13] Sam Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks, 2019.
[14] Johannes Hendriks, Carl Jidling, Adrian Wills, and Thomas Schön. Linearly constrained neural networks, 2020.
[15] Valerii lakovlev, Markus Heinonen, and Harri Lähdesmäki. Learning continuous-time pdes from sparse data with graph neural networks, 2020.
[16] Xiaowei Jia, Jared Willard, Anuj Karpatne, Jordan Read, Jacob Zwart, Michael Steinbach, and Vipin Kumar. Physics guided rmns for modeling dynamical systems: A case study in simulating lake temperature profiles. Proceedings of the 2019 SIAM International Conference URL http:
on Data Mining, pp. 558-566, May 2019. doi:
/ //xx. doi. org/ 10. 11111117711 1...77....1111. 978111191111111771
[17] Chiyu "Max" Jiang, Karthik Kashinath, Prabhat, and Philip Marcus. Enforcing physical constraints in (enn)s through differentiable (pde) layer. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2020. URL net/ / forum? QoUgUkKK
[18] Anuj Karpatne, William Watkins, Jordan Read, and Vipin Kumar. Physics-guided neural networks (pgnn): An application in lake temperature modeling, 2017.
[19] K. Kashinath, M. Mustafa, A. Albert, J-L Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzade- nesheli, R. Wang, A. Chattopadhyay, A. Singh, A. Manepalli, D. Chirila, R. Yu, R. Walters, B. White, H. Xiao, H. A. Tchelepi, P. Marcus, A. Anandkumar, P. Hassanzadeh, and null Prabhat. Physics-informed machine learning: case studies for weather and climate mod- elling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and En- gineering Sciences, (tfeences379(1194))))) 379(2194): 2021. doi: 12.1098/rtt.222200000 //
URL https: wiw.issiinggg..... doi/ abs / 10. 10w.1009///0 www.:/9//r//eeeee 2020.
[20] Y. Khoo, J. Lu, and Lexing Ying. Solving pde problems with uncertainty using neural-networke 2018.
[21] Byungsoo Kim, Vinicius C. Azevedo, Nils Thuerey, Theodore Kim, Markus Gross, and Barbara Solenthaler. Deep fluids: A generative network for parameterized fluid simulations. Computer Graphics Forum, 38(2)599-70, May 2019. ISSN 1467-8659. doi: 10. 10. 11111/ 13619. URL
http: / dx doi org/ 10. 1 111 111 /tt.////ffff./ 13619.
[22] Jiachen Li, Hengbo Ma, and Masayoshi Tomizuka. Conditional generative neural system for probabilistic trajectory prediction, 2019.
[23] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data, 2017.
[24] Gaurav Manek and J. Zico Kolter. Learning stable deep dynamics models, 2020.
[25] Arvind T. Mohan, Nicholas Lubbers, Daniel Livescu, and Michael Chertkov. Embedding hard physical constraints in neural network coarse-graining 3d turbulence, 2020.
[26] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style. high- performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alché- Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 8035- Curran Associates, Inc., 2019. URL 901 2015-pytorch-ane http://wwwc.cooeeeeeeeeeeeeeeone pdf.
[27] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez- Gonzzllzz and Peter W. Battaglia. Learning mesh-based simulation with graph networks, 2020.

---

Under review as a conference paper at ICLR 2022
[28] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Thysics-informet neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[29] Martin Riedmiller and Heinrich Braun. Rprop- fast adaptive learning algorithm. Technical report, Proc. of ISCIS VII,, Universitat, 1992.
[30] L. V. Rueden, Sebastian Mayer, Katharina Beckh, B. Georgiev, Sven Giesselbach, R. Heese, Birgit Kirsch, Julius Pfrommer, Annika Pick, Ramamurthy, Michal Walczak, Garcke, C. BBuuchhage, and Jannis Schuecker. Informed machine learning taxonomy and survey of integrating knowledge into learning systems. arXiv: Machine Learning. 2019.
[31] Alvaro Srrccezz-onzzlezz Nicolas Heess, Jost Tobias Springenberg, Josh Merel, Martin Ried- miller, Raia Hadsell, and Peter Battaglia. Graph networks learnable physies engines for inference and control. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th Inter- national Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 4470 4479. PMLR, 10-15 Jul 2018.
[32] W.E. Schiesser. The Numerical Method of Lines: Integration of Partial Differential Equations. Elsevier Science, 2012. ISBN 9780128015513, URL https: ptttts://boooe /tt.ks/oooos///oo/oooo.ooooosos. google..o.ooooooooooooooooo..... fi/
[33] Sungyong Seo and Yan Liu. Differentiable physics-informed graph networks, 2019.
[34] Naoya Takeishi and Yoshinobu Kawahara. Learning dynamics models with stable invariant sets, 2020.
[35] T. Tieleman and G. Hinton. Lecture 6.5 Rr--Propp Divide the gradient by running average of recent magnitude. COURSERA: Neural Networks for Machine Learning, 2012.
[36] Jared Willard, Xiaowei Jia, Shaoming Xu, M. Steinbach, and Kumar. Integrating physics- based modeling with machine learning: A survey. ArXiv, abs/2003.044993 2020.
[37] Zeng Yang, Jin-Long Wu, and Heng Xiao. Enforcing deterministic constraints generative adversarial networks for emulating physical systems, 2019.
[38] Haijun Yu, Xinyuan Tian, Weinan E, and Qianxiao Li. Onsagernet: Learning stable and interpretable dynamics using a generalized onsager principle, 2020.
12

---

Under review as a conference paper at ICLR 2022
A GENERALIZED CHAIN RULE AND HANDLING MIXED PARTIAL DERIVATIVES.
y g(1,. Fnn with all arguments being either identical, distinct or grouped. Then, partial derivatives of (V) can be evaluated using the Faà di Bruno's 's formula
8th f(v)
alBly
Σ TEM
(v)
II
0x Ozn
IIjEn dr,
BET
where 1I is the set of all partitions of the set {I....,,) B runs through elements of the partition f(() denotes m tt derivative, and | is cardinality.
n).
The formula consists of two terms: f(x (y), which can be evaluated using automatic differentiation, ø1 and I Smj" which can be evaluated as shown Equation 8.
In case that all F Zn are identical, the mixed derivative 8" f(x) reduces ("v(()) 0z
dzf
B DATA GENERATION
In all cases we run simulation fine grid and then interpolate the results to coarser grid represented as the "full grid" the experiments.
B.1 CAHN-HILLIARD EQUATION WITH AN INTEGRAL CONSTRAINT Training and testing data was obtained by solving
Ou
2V (u(1 - -)² -u) V
(14)
Öt
on unit interval with periodic boundary conditions and e = 0.04. The domain was represented by uniform grid with 100 nodes and the time step was 1.0e-6 sec. The initial conditions ug (z) were generated as follows
uo (x)
(A cos ((x 8)2m) + 7 sin ((x - 8)2m))
(15)
(x) min un (x)
un (x) max (z min uo (x))
(16)
where ds. Unif(-1,) and s Onif(0, 1). Examples of the simulations are shown in Figure 10.
B.2 HEAT EQUATION WITH MONOTONICITY CONSTRAINT
Training and testing data obtained by solving
Ou.
DV
(17)
Ôt
on a unit square with zero Neumann boundary conditions and = 0.2. The domain was represented by unstructured grid with 2971 nodes and the time step was set 0.001 sec. The initial conditions (x)) were generated as
[wx)
Σ war"
(18)
1 3
(), cos ((x + s)2m) + 7/4 sin ((x s)2x)) + Ao
((v) 1/2
Σ i=1
(19)
2
uo(x, 8) f(x) + g(v),
(20)
uo y) min (2 v)
uo(a, 4) max to 3) min on (x. 8) where w; ~ 1ni 1.1), Ao, 8 Unif(-1,1)
(21)
Examples of the simulations shown in Figure/ 11.
13

---

Under review as a conference paper at ICLR 2022
0 sec
0.0015 sec
0.0030 sec
Figure 10: Examples of simulations for the Hahn-Hilliard equation on a unit interval.
B.3 GAN WITH A DIVERGENCE CONSTRAINT
The data was generated by sampling random velocity fields and then projecting them to the space of divergence-free fields. The procedure was as follows. First, a random velocity field uo y) was generated on unit square by generating each component i
v) = Σ
Aua cos (ko + lyy + Ya sin (ker + Iy),
(22)
k,==-N
Lox( v) min (xx v)
uns (x, v) = 6 max ūox (o; y) min (o:( v)
0.5
(23)
where N 10 and λRa, Mat (0, 1). Then, the divergence-free component of uo (z, y). denoted
by un (z, y). was extracted by using the projection method by solving V Ug = V2 @ for and then evaluating unf (z, y) un (z, y) Vg. Finally, the data was scaled to [11,
B.4 GAN WITH A LAPLACIAN CONSTRAINT The data was generated by solving
V
(24)
on unit square with Dirichlet boundary conditions. The domain was represented by an unstructured grid with 2971 nodes. The boundary conditions were generated by generating random functions up (z) and using their boundary values the boundary conditions. The functions uo were generated as
uz 3) =
Σ
Ae cos (kx + Ig) + yuu sin (kr + ly)
(25)
k.1--N N
where N 5 and Aut, Tha ~ NO, 1). The data was then scaled [0, 1].
14

---

Under review as a conference paper at ICLR 2022
sec
0.1 sec
0.2 sec
Figure 11: Examples of simulations for the heat equation on unit square.
C MODELS, TRAINING AND TESTING
C.I REPLACING EXISTING METHODS
For our comparisons considered experiments from two works. Next, we provide some details about these experiments.
The first experiment was taken from [37] Section 3.2. The experiment shows how soft physics- constraints affect predictions of GAN learning distribution of divergence-five fields. The data generated on a uniform grid which allows to evaluate divergence using finite differences. The constraint is enforced through extra loss term which penalizes violation of the constraint. The performance metric used is the Frobenius norm the divergence averaged over all fields in batch. For training we used code provided by the authors with original parameters. We replaced finite differences in the constrained evaluation function with our method.
The second experiment was taken from [10]. This work deals with the task of predicting sea surface temperatures future times given snapshots of the temperature over current and previous times. The model proposed by the authors accomplishes this tasks by taking sequence of surface temperatures times fi.k..... 1; and predicting the underlying motion field which is then used to predict the temperature time f.+1. Insights about physical properties of the motion field were used to constrain the model's predictions. Constraints are imposed on divergence, magnitude and gradients of the motion field. The data is generated on uniform grid which allows evaluate the constraints using finite differences. The constraints are enforced through extra loss terms which penalize violation of the constraints. Performance metrics that were used are MSE between the data and model predictions, smoothness loss and divergence loss. For training we used code provided by the authors with original parameters. We replaced finite differences in the constrained evaluation function with method.
15

---

Under review as a conference paper at ICLR 2022
C.2 CAHN-HILLIARD EQUATION WITH AN INTEGRAL CONSTRAINT
In all experiments with the Cahn-Hilliard equation represent the dynamics function Fa by MLP with 2 hidden layers and LeakyReLU nonlinearities (negative slope 0.2). The number of neurons in each layer was set to the number of nodes in the spatial grid which the model was trained. The predictions u(t) were obtained by simulating the system forward in time using adaptive Heun solver from torchdiffeq package [6] with rtol and atol set 1.0e-5 and 1.0e-5 respectively. All models were trained for 1500 epochs using Rprop optimizer [29] with learning rate set to 1.0 10 and batch size set to the number of simulations in the training set. Mean squared error was used as the loss function. Spatial and temporal grids in the testing data were the same as in the training data. We set À = 2.
C.3 HEAT EQUATION WITH A MONOTONICITY CONSTRAINT
In experiments with the heat equation we represent the dynamics function Fe by a GNN with the messaging and aggregation networks being MLPs with single hidden layer consisting of 60 neurons with Tanh nonlinearities and the input/output sizes of 4/40 and 41/1 respectively. The predictions u(t) were obtained by simulating the system forward in time using adaptive Heun solver from torchdiffeq package [6] with rtol and atol set to 1 1.0e-5 and 1.0e-5 respectively. All models were trained for 750 epochs using Rprop optimizer [29] with learning rate set 1.0 10-6 and batch size the number of simulations in the training set. Mean squared error was used as the loss function. Spatial and temporal grids in the testing data were the same in the training data. We À 0.1.
C.4 LEARNING DISTRIBUTIONS OF PHYSICAL FIELDS
In both cases used identical architectures and training process for the constrained and uncon- strained models. Both models were trained for 1.2M iterations using the same random seed. Con- straints in the constrained model were enabled only after 600k iterations. The base distribution was set to a 228-dimensional isotropic standard normal. Models were trained using RMSProp optimizer [35] with batch size and learning rate set to 64 and 0.00001 respectively. The discriminator's weights were clipped to [-0.00 0.001, 0.01].
C.4.1 ZERO-DIVERGENCE FIELDS
We used MLPs as a discriminator and generator. The discriminator consisted of 3 hidden layers of sizes 1024-512- 256 with LeakyReLU nonlinearities (negative slope 0.2) and input/output size of 2010 and respectively. The generator consisted of 3 hidden layers of sizes 256-512-1024 with LeakyReLU nonlinearities (negative slope 0.2), input/output size of 128 and 2010 respectively, and final hyperbolic tangent nonlinearity applied to the output. We set X = 0.2.
C.4.2 ZERO-LAPLACIAN FIELDS
We used MLPs as a discriminator and generator. The discriminator consisted of 3 hidden layers of sizes 1024-512-256 with LeakyReLU nonlinearities (negative slope 0.2) and input/output size of 1086 and respectively. The generator consisted of 3 hidden layers of sizes 256-512-1024 with LeakyReLU nonlinearities (negative slope 0.2), input/output size of 128 and 1086 respectively, and sigmoid function applied to the output. We set À = 0.0075.
D CAHN-HILLIARD EQUATION WITH HARD INTEGRAL CONSTRAINTS
Here we demonstrate how the approach to enforeing hard constraints described in Section 3.3 can be used enforce integral constraints on a nonuniform grid.
We use the same setup as in Section 4.2 with 30 training simulations and 50% of nodes in the grid. We compare three models: unconstrained model, model with soft constraint and model with hard constraint. We use the PWL basis during training and testing.
Table shows that relative errors of all three models are practically similar but, as Figure 12 demonstrates, constraint violations differ significantly. We see that constraint violations of the model with hard constraint are zero at all time points as expected.
Being able to produce predictions that satisfy some constraints exactly might be very useful for some applications. however, constraints is limited to systems with a relatively small number of nodes and is significantly slower than models with soft constraints. We report training times in Table 1.
mention in Section 3.3, currently this approach to enforcing hard
16

---

Under review as a conference paper at ICLR 2022
Table
Test relative errors and training times for the
equation.
Model
Test rel. err
Training time
Unc.
0.041 + 0.002
15 min.
Soft-con.
0.044 + 0.002
17 min.
Hard-con.
0.043 + 0.002
207 min.
Unc.
Soft
Hard
0.010
Train
0.010
Train
0.010
Train
Conetr violation
Const: violation
Constr violation
0.005
0.005
0.005
0.000
0.000
0.000
0.000
Time (sec)
Time (sec)
Time (sec)
Figure 12: Constraint violation plots for the Cahn-Hilliard equation. Each panel shows the mean absolute violation of the constraint, Sa uy (1, CJ2 over all train/test simulations for each time point.
E LEARNING
EQUATION WITH GNNS
We the same setup in Section 4.2 with 75% of nodes in the grid but trained the models for 1500 epochs. Instead of an MLP we use a GNN with messaging and aggregation networks being MLPs with two hidden layers size 64 and LeakyReLU nonlinearities (negative slope 0.2). For each node, the GNN evaluates the output as
du
Σ JEN(S
po proj
(26)
where w!"T and mi are linear projections of the state at nodes and j and æ.y is a linear projection of the pair consisting of the distance between nodes and j and a unit vector pointing from node to node All projections have dimension 16.
compare constrained and unconstrained models. We use the PWL, PWQ and cubic RBF bases. Results of the experiment shown in Figure 13. The figure shows that relative errors and constraint violations of all models are significantly higher than for MLP-based models.
17

---

Under review as a conference paper at ICLR 2022
Node frac. 0.75
0.20
Rel, error
I
pwq
0.15
pwl
cubiC
0.10
30
60
120
0.050
Constr. viol.
0.025
0.000
30
60
120
Train size
Figure 13: Relative errors and constraint violations for a GNN trained on the Cahn-Hilliard equation. All resuts are for the test set.
EXTRA FIGURES
1.
02
Figure 14: Grid used for training GANs.
18


Reviews:
[{"review_id": "lwUITMyUF2m", "reviewer": "Reviewer_A9JW", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "In general, the paper is easy to read. Many details remain vague though, e.g., the implementation of the constraints in the experiments, how to train the differentiable hard constraint solver, how exactly are constraints split into gamma_n terms (and why, this seems unnecessary), what is a thin plate spline and why is this a good basis for evaluating constraints.\nSection 3.1, 3.2 could have been framed better by noting that this is simply a finite element formulation of constraints, and differential/integral formulations directly correspond to FE differentials/integrals.\n\nAs the paper mentions in the introduction, using soft constraints (often on regular grids) in training is an extremely common practice in physics prediction papers. Often, this is not referred to as constraints but as \"physics-based\"/auxiliary losses (e.g. adding a loss on divergence when modeling incompressible flow), but it amount to the exact same thing.\nThis paper differentiates itself from these methods by a) using unstructured grids, and b) providing a formalism for solving hard constraints.\n\nHowever, the paper doesn't make a particular strong case for these choices. First, the results are either on regular grids (4.1), or irregular meshing of a square domain, but with a very uniform edge lengths (4.2-4.3). None of these demonstrate any advantage of irregular meshing, and I expect all examples could have been much more easily solved on regular grids with standard CNN-based methods.\nAnd second, hard constraints are presented as one of the contributions, and the description spans quite a bit of the method section. However as the authors state themselves, their method for solving hard constraints is only feasible for very small systems, and is hence not used for any of the paper's main results (only for a minor example in the appendix).\nSimilarly, the paper introduces mitigation techniques to make quadratic basis functions (PWQ) work, which turn out to perform strictly worse than a simple linear basis in all settings. \n\nAnd finally, I'm a bit unsure what to take away from the results. It's a bit of an obvious finding that putting an auxiliary loss on a certain quantity (e.g. divergence) will produce results with lower values for this quantity. It would have been much more interesting to show examples that actually highlight the contributions of the paper-- i.e. that isn't easily achieved with corresponding loss terms on regular grids. Also, this paper really needs comparisons to other methods, e.g. PDEs on regular grids with physics-based losses, or perhaps comparing constraints on energy to Hamiltonian methods, etc.\n\nOther comments:\n- On regular grid in 4.1, wouldn't the used PWL basis be exactly equivalent to FD?\n- While there seems to be a bit less noise, I'm not sure that the constrained look 'more similar' to the data in fig 6, 9\n- Why does the GNN in 4.2 not work? GNNs generally work very well for such PDEs. An MLP is likely a poor choice for learning local PDE dynamics, particularly when trained on small datasets", "summary_of_the_paper": "The paper introduces and studies several variants of introducing soft and hard constraints into learned models of PDEs on unstructured meshes.", "summary_of_the_review": "While I like aspects of this paper (e.g. I think the FE formulation for constraints on unstructured grids can be useful), in its current form it's not ready for publication. Many of the more interesting contributions (e.g. higher-order basis functions, hard constraints, irregular grids) don't really work that well and aren't used in the main experiments: most results are on square, regular domains, with soft constraints and linear or radial basis functions, which is quite similar to what a lot of physics prediction papers are already doing with auxiliary losses (and for which there aren't any baseline comparisons). \nI'd encourage the authors to work further on the method (e.g. make hard constraints work for nontrivial systems), and choose good experiments to showcase the strengths of constraints on irregular grids. I also think constraints for adversarial methods that the paper touches on could be worth exploring further.", "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": {}}, {"review_id": "lYqtJ_tZqF", "reviewer": "Reviewer_vTXr", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "Advantages\n\n- The submission is written well, is clear up until the experimental section.\n\n- The experiments show that the addition of constraints help to learn a model that does not violate said constraints.\n\n- The central method is sound and simple enough to be reimplemented. That said, I think some experimental details are missing and I could not reimplement the experiments. That said, since code will be released, this is not such a big issue.\n\nQueries\n\n- Equation 1: can dynamics F also explicitly depend on t? As written, it does not.\n\n- Equation 7: I’m confused about the dimensionality of \\alpha_j. In Equation 6, it is a p-dimensional quantity, because u(t, x) \\in R^p. Is it the case that u(t, x) is now scalar-valued, since h: R \\to R and \\gamma_m: R \\to R? This should be made explicit in the text.\n\n- Equation 10: what is g?\n\n- Section 4.1: I think it would be more helpful to the reader if the experimental setups from [37] and [10] were described more completely in the main text. It is difficult to gauge from what is written, what exactly is going on and whether your results are meaningful\n\n- Section 4.2: Could you provide intuition for why the MLP would work better than the GNN? I would hazard a guess that the GNN does not have a large enough receptive field. Since you only use a single layer of message passing from what I see in the appendix, the receptive field of an output neuron is the 1-step neighborhood due the the Delaunay triangulation. By contrast, the MLP receives the full domain as its receptive field.\n\n- There seems to be a lack of baselines in the experimental section after section 4.1. I am not entirely sure why this is. Given that it is mentioned in the related work that there is a rich literature on incorporating constraints, both hard and soft, into learnable PDE models, I would expect to have seen some in the experiments.\n\n- What do the error bars represent in Figures 4 and 7?\n\nMinor notes\n\n- Equation 2: the ordering of x and t in du(x, t)/dt does not align with elsewhere in the submission, where you have used u(t, x).\n\n- Equation 2: Please be more specific what the notation F_\\theta(x_i, x_N(i), …) means. For instance what is x_N(i) or u_N(i)? My understanding is that you have replaced the differential operator F with a local parametrised operator F_\\theta over neighbourhoods N(i). The notation and underlying assumptions here should be explained more precisely.", "summary_of_the_paper": "The authors present a method to incorporate constraints into the output of learnable PDE models. They cover point-wise, differential, and integral constraints. They achieve this by representing the PDE solution in a basis, as is common for the variational method (e.g. pseudo-spectral method, and finite element method). The neural network outputs interpolation coefficients for each time step into the future. To implement the constraints, they have two methods: ‘soft constraints’, whereby a constraint-breaking penalty is applied to the neural network output and ‘hard constraints’, whereby the output of the network is projected on to constraint-satisfying solutions. This latter method is achieved by solving a convex programme.\n\nThe authors test on a variety of benchmarks, demonstrating that their method works; although, it is hard to parse the experimental section for whether these results are significant.", "summary_of_the_review": "I think the method is sound, and up until the experimental section the paper is well laid out and straight-forward to read. My understanding began to breakdown in the experimental section, where information was either referenced in the Appendix or completely in other papers. This made it hard to read and slowly down reviewing a lot. In most of the experiments, comparisons to other methods in the literature are missing and it is therefore hard to gauge the significance of this work.", "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": {}}, {"review_id": "EAE5BXl-7BX", "reviewer": "Reviewer_3bQB", "summary": "", "questions": "", "limitations": "", "rating": 5, "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 task addressed by the authors is crucial in machine learning since purely data driven approaches blatantly fail at conserving several physical quantities such as momentum or energy. In particular, the enforcement of algebraic constraints (whether differential or integral) on a learned model is a necessary path to learn from/for physical data.\n\nThe method can be summarized as follows: Learn the model $F_\\theta$ that fits the data constraints (regression, gan etc. ...), hence the solution at “grid points” $\\mathbf{u}$. Then fit a simple function of separate space and time $u_f (t,x) = \\alpha_j(t) \\phi_j(t)$. (linear, quadratic, etc..) to interpolate between the prediction on which can enforce specific constraint, then a fortiori constraining the learned $u$. The experiments seem to support the author’s choices. The experiments are well conducted and present various use cases useful for the whole community.\n \nHowever, the main weakness of the paper is the clarity in the presented interpolation method. Mainly, how it is computed and its link to the learned solution $\\mathbf{u}(t)$. Also the estimation of $\\phi$ is only very lightly discussed.\n\nQuestion:\n\n1. How well can the interpolant $u_f$ approximate the original to the problem ? Indeed, I guess this is related to the convergence of a numerical scheme: if the space is well covered (i.e. discrete step size goes down to 0) can we expect to recover a good function $u_f$ using such a prior?\n2. Why piecewise cubic approximations do not require extra conditions (do you use Hermite polynomials and continuity of the derivative) compared to piecewise quadratic ? I suggest the authors include a thorough discussion on this topic in the appendices (with the conditions imposed on such local polynomials). \n3. For the gan experiments, does the generator depend on time ? If  not, then the interpolation simply concerns the $\\phi$ part ?\n4. How to enforce a second order differential constraint when the interpolation is PWL ?\n\nRemark on the presentation:\nDespite interesting ideas, there is room for improving the presentation and the clarity of the paper.\n\na) Since the main focus of the experiments concerns “soft constraints”, I suggest the authors develop the section treating the soft-constraint components (Perhaps including an algorithm to show practical use) and how constraining $u_f$ constrains the learned $\\mathbf{u}$ (or more specifically how the losses on $u_f$ impacts the $\\theta$ of $F_\\theta$), by detailing explicitly the link between the interpolation and the learned solution.\n\n\nb) The key point of the paper lies in the interpolation scheme, however, in my opinion the paper lacks elements of details for an in depth understanding of the estimation method of the parameters $\\alpha_j$ and $\\phi$. Since this point is central in the paper, I highly recommend including thorough details on the estimation of the $\\alpha$ and $\\Phi$ at least in the appendices.\n\nMinor:\n - I suggest for clarity to write the learned solution $u_\\theta$ instead of $u$ to avoid the possible confusion between $u$ the solution to (1) and $\\mathbf{u}$ the learned numerical solution at the nodes $x_1, .., x_N$). \n- Should $t$ be a variable of $F_\\theta$ in eq.2. Does $u_1 = u(x_1)$. ?\n- Eq.6. If there are $N$ points ($x_1$, …, $x_n$). Should the sum stop at $N-1$ ?\n- P.6 “with a multiplayer perceptron”\n", "summary_of_the_paper": "The paper proposes a two-folded method to enforce constraints of different natures (differential, integrals….) on a statistical model learned from physical data.\nThe constraints are not enforced directly on the model but rather on “interpolant” functions that aims at completing the original model in between the observed grid points.\n", "summary_of_the_review": "The paper proposes a very interesting method to enforce differential / Integral and differential constraints on an interpolation function derived from the learn operator that fits the data. \nThe experiments are well conducted but the method’s presentation lacks clarity (see main review) notably on the estimation of the interpolation coefficients and their link to the learned solution function? \n\nIf my concerns are addressed, I will be glad to increase my score.\n", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "9-rItcnq9zo", "reviewer": "Reviewer_HR2b", "summary": "", "questions": "", "limitations": "", "rating": 5, "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 presents a series of tests which follow a somewhat unintuitive order. First, a regular grid case is presented, which seems to be taken from previous work. Not much detail is given, and the test only confirms that the constraints give very similar results. This test seems like a pure debugging case, as it is based on Cartesian grids instead of the unstructured grids which the paper wants to focus on.\n\nThe subsequent sections (4.1 ff) then provide examples with unstructured grids. While the first CH-case is very simplistic, the latter cases contain interesting setups and show some interesting results. E.g., I found it interesting to see how the physical constraints affect the smoothness of the zero divergence fields \n\nIn terms of writing, I found the order and argumentation of the results sub-optimal, as mentioned above. The structured grid case seems mostly out of place. The introduction also sounds strange to me, with current work being \"limited\" to uniform grids. Previous works have admittedly focused on these, but the methods are still applicable - after all, this submission is applying the methods from uniform grids in a fairly straight forward manner. I would recommend to rephrase this. Also, stylistically, citations shouldn't be used as nouns, and the related work seems to brief to me (e.g., generative models for physical problems are missing).\n", "summary_of_the_paper": "This paper proposes a method to enforce physical constraints in deep learning models. It provides a nice summary of local, differential and integral constraints, and frames them in a Lagrangian setting. For a reason which I could not clearly follow, the paper focuses on GAN early on. This is intuitive to me, as the physical constraints could nicely stand on their own. \n", "summary_of_the_review": "Overall, I think this submission would benefit from further revisions. It is definintely an interesting direction, but in its current state the paper applies existing methods from regular grids to Lagrangian discretizations. This is a good idea, but not a very fundamental step forward. In combination, the results do not contain many surprises, and the focus on GANs in a way detracts from the main goals of the paper. The presentation of the work could also be further improved. So overall, I'm leaning towards the negative side given the current state of the submission.\n", "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 introduces a framework for enforcing constraints into deep NNs used for modeling spatio-temporal dynamics characterizing physical systems. The authors consider different types of constraints (pointwise, differential and integral). They start from a formulation approximating PDEs as set of ODEs (method of lines). Their main idea is to approximate the solution of the equations using an interpolant between observations and imposing the constraints on this approximation function. The interpolant is built using basis functions located at observation points. The formalism considers irregular spatial grids and both soft and hard constraints. The main claim is then the introduction of a general formalism for considering different types of constraints on irregular grids. Experiments illustrate the behavior of the proposed method on different types of evolution equations and constraints.

The reviewers agree that the proposed approach is interesting and that some of the ideas are original. However, they also consider that the paper is not convincing enough to demonstrate the interest and novelty of the approach, compared to alternative methods. The experimental section mainly considers (except for one application) regular grids and constraints that could be handled by other methods as well. The authors should present cases where their method provides a clear advantage, distinct from existing solutions. The authors provided a well-argued rebuttal, clarifying several points. However, all reviewers retained their original scores and encourage the authors to further develop the experimental analysis to present a stronger paper. In addition, the presentation could be improved, and some technical aspects better explained (e.g., description of interpolation methods, and some advice on which interpolant to choose for a given problem).
Decision: Reject