Year: 2021

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 a model to satisfy a set of equations representing physical constraints. Currently available approaches can enforce a very limited set of constraints and are applicable only to uniform spatial grids. We propose a method for enforcing general pointwise, differential and integral constraints on unstructured spatial grids. Our method is based on representing a model’s output in terms of a function approximation and enforcing constraints on that approximation. We demonstrate wide applicability and strong performance of our approach in data-driven learning of dynamical PDE systems and distributions of physical fields.

heading: 1 INTRODUCTION
text: Multiple works have shown the capability of neural networks to 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 physics-constrained 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 we consider in this work is based on forcing a model to satisfy algebraic constraints represented by a set of equalities and inequalities. This is the most commonly used approach which allows to represent a wide range of constraints and has been shown to work well in many cases [18; 17; 25]. However, while many constraints can be represented algebraically, it is not always clear how to evaluate and enforce them.
Currently available approaches to enforcing algebraic constraints are 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 can be readily applied to models based on convolutional neural networks (CNNs) but cannot be extended to recently developed models based on graph neural networks (GNNs) [33; 27; 15] and other models working on unstructured grids.
We propose a much more general method which allows to 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 a grid to construct an interpolant and applying constraints directly to that interpolant (Section 3).
Code and data will be made publicly available.

heading: 2 BACKGROUND
text: PDE-driven dynamical systems. Many physical systems can be described in terms of PDEs. Such systems are defined on a bounded domain on which they evolve over time. We consider continuous dynamical systems with state u(t,x) ∈ Rp that evolves over time t ∈ R≥0 and spatial locations x ∈ Ω ⊂ 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 an unknown PDE
∂u(t,x)
∂t = F (x, u(t,x),∇xu(t,x),∇2xu(t,x), ...) (1)
which describes the temporal evolution of the system in terms of the locations x, state u and its first and higher-order partial derivatives w.r.t. x. The goal of a 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 (x1, . . . ,xN ) over increasing time points (t0, . . . , tM ). This results in a dataset (y(t0), . . . ,y(tM )), where y(ti) = (u(ti,x1), . . . , u(ti,xN )) is a collection of observations. The dataset is used to train the model to predict (y(t1), . . . ,y(tM )) starting from the initial state y(t0). Training is typically done by minimizing an average loss between the model’s predictions u(t) and the data y(t).
PDE models differ in 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 1 which results into a system of ODEs
u̇(t) :=  du(t,x1) dt ...
du(t,xN ) dt
 ≈  Fθ(x1,xN (1), u1, uN (1))... Fθ(xN ,xN (N), uN , uN (N))  (2) which approximates the solution of Equation 1 at the observation locations xi using their neighboring points N (i), where xN (i) and uN (i) are the neighbors’ positions and states respectively, and ui is u(t,xi). The approximate solution converges to the true solution as N increases. The true dynamics F is approximated by a parametric model Fθ whose parameters θ are learned by minimizing the difference between the model’s predictions
u(t) = u(0) + ∫ t 0 u̇(τ)dτ (3)
and the data y(t). The integral in Equation 3 is solved using a numerical ODE solver. In [15], the function Fθ was represented by a graph neural network (GNN) which takes states and locations at an observation point i and its neighboring points N (i). The observation points are connected into a grid using Delaunay triangulation which allows to naturally defineN (i) as a 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 θ are learned by minimizing the MSE between y(t) and u(t)
Ldata = 1
M M∑ i=1 ‖u(ti)− y(ti)‖22. (4)
The gradient of Ldata w.r.t. θ is evaluated using the adjoint method as 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 a generative model consisting of a generator and a discriminator [12]. The generator, G, learns to transform a random variable Z ∼ pZ over a latent space Z to the data space Y in such a way that the discriminator, D, cannot tell the difference between samples generated by G and samples from the data distribution pdata. Both, G and D are learned by solving the following minimax problem
min G max D V (G,D) = EY∼pdata [logD(Y )] + EZ∼pZ [log (1−D(G(Z)))] . (5)
Solution of this problem exists and is unique with the optimal generator perfectly mimicking the data distribution [12].

heading: 3 METHODS
text: 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.

heading: 3.1 EVALUATING CONSTRAINTS ON UNSTRUCTURED GRIDS
text: We assume the data y(t) is available at observation points (x1, . . . ,xN ) and time points (t1, . . . , tM ) and that a 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 on it directly. Instead, we approximate it by an interpolant uf (t,x) and impose constraints on uf (t,x) (Figure 1). The approximation is constructed from u(t) by placing a basis function at each xi and representing uf (t,x) as
uf (t,x) = N∑ j=1 αj(t)φj(x), (6)
where φj is a scalar basis function at xj and αj ∈ Rp. The coefficients αj(t) are obtained from u(t) (see Section 3.4).
Next, we show how to evaluate constraints on uf (t,x) using basic building blocks. To avoid cluttered notation, we consider equality constraints and assume u(t,x),x ∈ R. Generalization to inequality constraints, vector fields and higher spatial dimensions is straightforward.
Pointwise constraints. Consider points z = (z1, . . . , zK) in Ω on which a pointwise constraint h(uf (t, zi)) = 0 should be evaluated. Assume the function h : R→ R is representable in terms of a finite number of functions γm(uf (t, zi)) : R → R indexed by m. For example, should the constraint be h(uf ) = 3uf + u2f = 0, then we
would define γ1(uf ) = uf , γ2(uf ) = u2f and h(uf ) = 3 · γ1(uf ) + γ2(uf ) = 0. Then, h can be evaluated by evaluating each γm as
γm(uf (t, zi)) = γm  N∑ j=1 αj(t)φj(zi)  = γm (Φi,·α(t)) , (7) where α(t) = (α1(t), . . . , αN (t))T , Φ is K-by-N matrix with elements Φi,j = φj(zi), and Φi,· is the i’th row of Φ.
Differential constraints. Consider the same setup as before but now h(uf (t, zi)) = 0 consists of differential operators and is representable in terms of a finite number of functions ∂ qγm(uf (t,zi))
∂zqi :
R→ R indexed bym, where the derivative order q could be different for eachm. For example, should the constraint be h(uf ) = 3uf + uf · ∂u2f ∂x = 0, then we would define γ1(uf ) = uf , γ2(uf ) = u 2 f and h(uf ) = 3 · γ1(uf ) + γ1(uf ) · ∂γ2(uf )∂z = 0. Then, h can be evaluated by evaluating each ∂qγm(uf (t,zi))
∂zqi using the generalization of the chain rule (Appendix A) which contains only two types
of terms. The first type of terms dγmduf , . . . , dqγm duqf can be evaluated using automatic differentiation while the second type of terms ∂uf∂zi , . . . , ∂quf ∂zqi can be evaluated as
∂quf ∂zqi = N∑ j=1 αj(t) ∂qφj(zi) ∂zqi = Φ (q) i,· α(t), (8)
where Φ(q)i,j = ∂qφj(zi) ∂zqi . Mixed partial derivatives can be handled in a similar way (Appendix A).
Integral constraints. Consider the same setup as before but with h(uf (t,x)) =∫ Ω τ(uf (t,x))dx = 0, where the function τ : R → R is representable in terms of functions
γm(uf (t, zi)) : R→ R similarly to the pointwise constraints. Then, ∫
Ω τ(uf (t,x))dx can be evalu-
ated using a numerical integration technique, e.g. midpoint rule, Gaussian quadrature or Monte-Carlo integration, as ∫
Ω τ(uf (t,x))dx ≈ K∑ i=1 τ(uf (t, zi))µi, (9)
where K is the number of integration points, µi are integration coefficients which depend on the grid and integration method, and τ(uf (t, zi)) is evaluated as in Equation 7.

heading: 3.2 SOFT CONSTRAINTS
text: Soft constraints are implemented by minimizing the following loss Ldata + λr(h(uf )), where λ ∈ R and Ldata is defined as in Equation 4. We set r(h(uf )) = 1KM ∑K i=1 ∑M j=1 h(uf (tj , zi)) 2 for point-
wise and differential constraints and r(h(uf )) = 1M ∑M j=1 h(uf (tj ,x)) 2 for integral constraints.

heading: 3.3 HARD CONSTRAINTS
text: Our method allows to implement hard constraints by projecting the interpolant uf (t,x) to a subset of functions which satisfy the required constraints. Namely, if uf (t,x) does not satisfy constraints g and h, it is projected to a subset of functions which satisfy the constraint by solving the following optimization problem
min ûf∈Vφ
‖uf − ûf‖2L2
s.t. h(ûf ) = 0, g(ûf ) ≤ 0,
(10)
where the projection is denoted by ûf (t,x) and Vφ is spanned by the basis functions. Using the basis representation uf (t,x) = ∑N i=1 αi(t)φi(x) and ûf (t,x) = ∑N i=1 βi(t)φi(x) we can rewrite the optimization problem (10) as
min β(t)∈RN
(α(t)− β(t))T Φ̂(α(t)− β(t))
s.t. h(ûf ) = 0, g(ûf ) ≤ 0,
(11)
where β(t) = (β1(t), . . . , βN (t))T and Φ̂i,j = ∫
Ω φi(x)φj(x)dx.
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 do that, we restrict ourselves to constraints that can be expressed as an equality or inequality between Aβ(t) and b, where A is a constant matrix and b is a constant vector. This formulation admits pointwise, differential and integral constraints on untransformed uf . The objective function is convex since its Hessian is positive-semidefinite i.e. for any v ∈ RN
vT Φ̂v = N∑ i,j=1 vivjΦ̂i,j = N∑ i,j=1 〈viφi, vjφj〉L2 = 〈 N∑ i=1 viφi, N∑ j=1 vjφj〉L2 ≥ 0. (12)
This allows to solve the problem (11) and differentiate its solution β∗(t) w.r.t. α(t). The model parameters are found by minimizing the following loss function Ldata +λLproj, where λ ∈ R and Ldata is defined as in Equation 4 but with u(ti) replaced by û(ti) = (ûf (ti,x1), . . . , ûf (ti,xN )). We set Lproj = 1NM ∑N i=1 ∑M j=1 ‖uf (tj ,xi)− ûf (tj ,xi)‖ 2 2. The second term makes the optimization procedure prefer models that predict uf close to the feasible set of the problem (11).
We note that the proposed approach is currently limited to small-scale problems due to existing computational bottlenecks in the implementation of differentiable convex optimization [1].
3.4 BASIS FUNCTIONS
Selecting appropriate basis is crucial for efficiency and applicability of the proposed method. Ideally, the basis should allow efficient construction of uf (t,x) from u(t), contain no tunable parameters, and lead to sparse matrices Φ,Φ(q) and Φ̂. We consider 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 to sparse Φ,Φ(q) and Φ̂. For the piecewise linear basis the interpolant uf (t,x) can be constructed directly from the predictions by setting α(t) = u(t). However, constructing uf (t,x) for a higher order basis, e.g. piecewise
quadratic, requires the model to make predictions not only at 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 α(t) 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 a wider range of properties. Some RBFs have tunable parameters, some don’t. The matrices Φ,Φ(q) and Φ̂ evaluated with RBFs are typically dense, but RBFs with compact support exist (e.g. bump function). The interpolant uf (t,x) can be constructed by evaluating α(t) = K−1u(t), where K−1 is the inverse of the interpolation matrix of the given RBF and Kij = φ(‖xi − xj‖), where φ is an RBF and xi,xj are observation locations. In this work we use the cubic RBF basis i.e. φ(r) = r3.
We use PyTorch [26] to handle sparse matrices and to evaluateK−1u(t) in a differentiable way.

heading: 4 EXPERIMENTS
text: In the following experiments we use the relative error between the data y(t) and model predictions u(t) defined as ‖y(t)−u(t)‖2‖y(t)‖2 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 a single NVIDIA Quadro P5000 GPU. All errors bars represent one standard deviation of the results over five random seeds.

heading: 4.1 REPLACING EXISTING METHODS
text: 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 of sea surface temperature by observing snapshots of the temperature over time while enforcing gradient and divergence constraints on the warping 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.
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 our approach produces a 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.

heading: 4.2 CAHN-HILLIARD EQUATION WITH AN INTEGRAL CONSTRAINT
text: We start with the 1D Cahn-Hilliard equation
∂u ∂t = 2∇2(u(1− u)2 − u2(1− u)− 2∇2u) (13)
which is known to conserve the state u i.e. h(u) = ∫
Ω u(t, x)dx− C = 0 at all time points, where C = ∫
Ω u(0, x)dx is a constant. This is an example of a 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 in the previous section with a single γ1 being the identity function. We use PWL, PWQ and cubic RBF bases and compare the results to an 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.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
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 a GNN produced poor results. For that reason we represented the function Fθ with a multiplayer perceptron (MLP) which updates the state of each node based on the states of all other nodes in the grid (results for a 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.
The coefficients α(t) 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 uf (t,x) using the thin plate spline basis (φ(r) = r2 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 to bases used for training.
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 PWL 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 is 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 to grids with a large number of nodes due to sparsity of the constraint matrices and efficient evaluation of α.
4.3 HEAT EQUATION WITH A MONOTONICITY CONSTRAINT
We impose constraints on a 2D system governed by the heat equation ∂u∂t = ∇2u 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 ∂u∂x ≥ 0. The constraint is evaluated as shown in the previous section with γ1 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
is 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 × 51 grid placed on top of the original grid.
To learn the dynamics of the system we use the model from [15] directly with the messaging and aggregation networks being MLPs with a single hidden layer consisting of 60 neurons with Tanh nonlinearities and the input/output sizes of 4/40 and 41/1 respectively.
During testing, we use predictions of the models to construct an interpolant uf (t,x) using the thin plate spline basis and evaluate the constraint on that interpolant. This allows to make a fair comparison between the unconstrained model and different bases.
Figure 7 shows results of the experiment. We observe that changing the node fraction equally increases relative errors of all models and has noticeable effect on constraint violations, especially for the unconstrained model. Constrained models tend to show slightly higher or 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 a 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.
Figure 6 shows qualitative difference between predictions of constrained and unconstrained models. It can be noted that predictions from the constrained model have noticeably smoother contours thus making the field more monotonous in the horizontal direction.

heading: 4.4 LEARNING DISTRIBUTIONS OF PHYSICAL FIELDS
text: We demonstrate the effect of adding constraints to a GAN when learning distributions of physical fields on unstructured grids. We use Wasserstein 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. Constraints are enabled only after 600k iterations. Constrained models are trained similarly to the unconstrained ones but with a modified generator loss defined as LG + λ ln (1 + LC), where LG is the standard generator loss and LC is the constraint-based loss. We defineLC as the mean value of h(u)2, where h is a constraint evaluated at the centroid of each cell in the grid.

heading: 4.4.1 ZERO-DIVERGENCE FIELDS
text: Divergence-free vector fields are often encountered in solutions of fluid dynamics problems. The divergence-free constraint on a vector field u(x, y) = (u1(x, y), u2(x, y))
T is defined as h(u) = ∂u1∂x + ∂u2 ∂y = 0. The constraint is enforced using the PWL basis. We generated a dataset with 10k divergence-free fields on an unstructured grid with 1050 nodes (Figure 14)
and used a WGAN to learn a 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 a 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.

heading: 4.4.2 ZERO-LAPLACIAN FIELDS
text: Fields with zero Laplacian represent solutions to some PDEs, for example the steady-state heat equation. The zero-Laplacian constraint on a scalar field u(x, y) is defined as h(u) = ∂
2u ∂x2 + ∂2u ∂y2 = 0.
The constraint is enforced using the cubic basis as the PWL basis has zero second derivatives
everywhere. We generated a dataset with 10k Laplacian-free fields on an unstructured grid with 1050 nodes (Figure 14) and used a WGAN to learn a distribution over such fields. Note that the generated fields are not entirely Laplacian-free due to discretization errors.
Results of the experiment are shown in Figures 9b and 8. Similarly to 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 8 shows that the constrained model generates fields that have much lower constraint violation and Laplacian distribution very similar to that of the data.

heading: 5 RELATED WORK
text: Soft constraints. Soft constraints are widely used due to being relatively easy to 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 to implementing hard constraints are 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; 9; 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 on a uniform grid might be limiting.
Constrained GANs Works such as [37; 17] showed how physics-based constraints benefit training and quality of the generated samples but are also limited to uniform grids.

heading: 6 CONCLUSION
text: We presented a general approach to enforcing algebraic constraints on unstructured grids and showed how it can be used to enforce 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 on various tasks.

heading: REPRODUCIBILITY STATEMENT
text: 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.

heading: A GENERALIZED CHAIN RULE AND HANDLING MIXED PARTIAL DERIVATIVES.
text: Let y = g(x1, . . . , xn) with all arguments being either identical, distinct or grouped. Then, partial derivatives of f(y) can be evaluated using the Faà di Bruno’s formula
∂nf(y) ∂x1 · · · ∂xn = ∑ π∈Π f (|π|)(y) ∏ B∈π ∂|B|y∏ j∈B ∂xj ,
where Π is the set of all partitions of the set {1, . . . , n}, B runs through elements of the partition π, f (m) denotes m’th derivative, and | · | is cardinality.
The formula consists of two terms: f (|π|)(y), which can be evaluated using automatic differentiation, and ∂
|B|y∏ j∈B ∂xj , which can be evaluated as shown in Equation 8.
In case that all x1, . . . , xn are identical, the mixed derivative ∂nf(y) ∂x1···∂xn reduces to ∂nf(y) ∂xn1 .

heading: B DATA GENERATION
text: In all cases we run simulation on a fine grid and then interpolate the results to a coarser grid represented as the "full grid" in the experiments.

heading: B.1 CAHN-HILLIARD EQUATION WITH AN INTEGRAL CONSTRAINT
text: Training and testing data was obtained by solving
∂u ∂t = 2∇2(u(1− u)2 − u2(1− u)− 2∇2u) (14)
on a unit interval with periodic boundary conditions and = 0.04. The domain was represented by a uniform grid with 100 nodes and the time step was set to 1.0e-6 sec. The initial conditions u0(x) were generated as follows
ũ0(x) = 10∑ i=1 (λi cos ((x− s)2π) + γi sin ((x− s)2π)) + λ0 2 , (15)
u0(x) = ũ0(x)−min ũ0(x)
max ũ0(x)−min ũ0(x) , (16)
where λi, γi ∼ Unif(−1, 1) and s ∼ Unif(0, 1). Examples of the simulations are shown in Figure 10.

heading: B.2 HEAT EQUATION WITH A MONOTONICITY CONSTRAINT
text: Training and testing data was obtained by solving
∂u ∂t = D∇2u (17)
on a unit square with zero Neumann boundary conditions and D = 0.2. The domain was represented by an unstructured grid with 2971 nodes and the time step was set to 0.001 sec. The initial conditions u0(x) were generated as
f(x) = 6∑ i=0 ωix i, (18)
g(y) = 1
2 3∑ i=1 (λi cos ((x+ s)2π) + γi sin ((x+ s)2π)) + λ0 2 , (19)
ũ0(x, y) = f(x) + g(y), (20)
u0(x, y) = ũ0(x, y)−min ũ0(x, y)
max ũ0(x, y)−min ũ0(x, y) , (21)
where ωi ∼ Unif(0.1, 1.1), λi, γi, s ∼ Unif(−1, 1). Examples of the simulations are shown in Figure 11.

heading: B.3 GAN WITH A DIVERGENCE CONSTRAINT
text: 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 u0(x, y) was generated on a unit square by generating each component i as
ũ0i(x, y) = N∑ k,l=−N λkl cos (kx+ ly) + γkl sin (kx+ ly), (22)
u0i(x, y) = 6× (
ũ0i(x, y)−min ũ0i(x, y) max ũ0i(x, y)−min ũ0i(x, y)
− 0.5 ) , (23)
where N = 10 and λkl, γkl ∼ N (0, 1). Then, the divergence-free component of u0(x, y), denoted by u∗0(x, y), was extracted by using the projection method by solving∇ · u0 = ∇2φ for φ and then evaluating u∗0(x, y) = u0(x, y)−∇φ. Finally, the data was scaled to [−1, 1].

heading: B.4 GAN WITH A LAPLACIAN CONSTRAINT
text: The data was generated by solving
∇2u = 0 (24) on a 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 u0(x) and using their boundary values as the boundary conditions. The functions u0(x) were generated as
u0(x, y) = N∑ k,l=−N λkl cos (kx+ ly) + γkl sin (kx+ ly) (25)
where N = 5 and λkl, γkl ∼ N (0, 1). The data was then scaled to [0, 1].

heading: C MODELS, TRAINING AND TESTING

heading: C.1 REPLACING EXISTING METHODS
text: For our comparisons we 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-based constraints affect predictions of a GAN learning a distribution of divergence-free fields. The data is generated on a uniform grid which allows to evaluate divergence using finite differences. The constraint is enforced through an extra loss term which penalizes violation of the constraint. The performance metric used is the Frobenius norm of the divergence averaged over all fields in a 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 at future times given snapshots of the temperature over current and previous times. The model proposed by the authors accomplishes this tasks by taking a sequence of surface temperatures at times ti−k, . . . , ti and predicting the underlying motion field which is then used to predict the temperature at time ti+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 a uniform grid which allows to 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 our method.

heading: C.2 CAHN-HILLIARD EQUATION WITH AN INTEGRAL CONSTRAINT
text: In all experiments with the Cahn-Hilliard equation we represent the dynamics function Fθ by an 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 on 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 to 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−6 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.

heading: C.3 HEAT EQUATION WITH A MONOTONICITY CONSTRAINT
text: In all experiments with the heat equation we represent the dynamics function Fθ by a GNN with the messaging and aggregation networks being MLPs with a 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.0e-5 and 1.0e-5 respectively. All models were trained for 750 epochs using Rprop optimizer [29] with learning rate set to 1.0 · 10−6 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 λ = 0.1.

heading: C.4 LEARNING DISTRIBUTIONS OF PHYSICAL FIELDS
text: In both cases we used identical architectures and training process for the constrained and unconstrained models. Both models were trained for 1.2M iterations using the same random seed. Constraints in the constrained model were enabled only after 600k iterations. The base distribution was set to a 128-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.01, 0.01].

heading: C.4.1 ZERO-DIVERGENCE FIELDS
text: 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 1 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 a final hyperbolic tangent nonlinearity applied to the output. We set λ = 0.2.

heading: C.4.2 ZERO-LAPLACIAN FIELDS
text: 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 1 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.

heading: D CAHN-HILLIARD EQUATION WITH HARD INTEGRAL CONSTRAINTS
text: Here we demonstrate how the approach to enforcing hard constraints described in Section 3.3 can be used to 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 1 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, as we mention in Section 3.3, currently this approach to enforcing hard 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.

heading: E LEARNING CAHN-HILLIARD EQUATION WITH GNNS
text: We use the same setup as 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 of size 64 and LeakyReLU nonlinearities (negative slope 0.2). For each node, the GNN evaluates the output as
dui dt = γ  1 |N (i)| N∑ j∈N (i) φ(uproji , u proj j , x proj ij ), u proj i  , (26) where uproji and u proj j are linear projections of the state at nodes i and j and x proj ij is a linear projection of the pair consisting of the distance between nodes i and j and a unit vector pointing from node i to node j. All projections have dimension 16.
We compare constrained and unconstrained models. We use the PWL, PWQ and cubic RBF bases. Results of the experiment are shown in Figure 13. The figure shows that relative errors and constraint violations of all models are significantly higher than for MLP-based models.

heading: F EXTRA FIGURES

Review summaries:
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.  The 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.

The paper proposes a two-folded method to enforce constraints of different natures (differential, integrals….) on a statistical model learned from physical data. The 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. 

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.