{'fraction_non_alphanumeric': 0.028193450444589027, 'fraction_numerical': 0.009759271307742356, 'mean_word_length': 4.865380538477846, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermodynamics-informed graph neural networks

In this paper we present a deep learning method to predict the temporal evolution of dissipative dynamic systems. We propose using both geometric and thermodynamic inductive biases to improve accuracy and generalization of the resulting integration scheme. The first is achieved with Graph Neural Networks, which induces a non-Euclidean geometrical prior with permutation invariant node and edge update functions. The second bias is forced by learning the GENERIC structure of the problem, an extension of the Hamiltonian formalism, to model more general non-conservative dynamics. Several examples are provided in both Eulerian and Lagrangian description in the context of fluid and solid mechanics respectively, achieving relative mean errors of less than 3% in all the tested examples. Two ablation studies are provided based on recent works in both physics-informed and geometric deep learning.


Introduction
Computer simulation has become a standard tool in a wide variety of fields, from natural to social sciences, in order to simulate physical phenomena and predict its future behaviour. Simulation models enable engineers to find optimal designs which can later be validated in a much more expensive experimental set-ups. Thus, it is important to develop models based on mathematical and physical foundations which ensure the reliability of the results. With the irruption of the information era, recent research lines have focused on the use of data-based deep learning algorithms which overcome some of the limitations of traditional methods, such as handling highly nonlinear dynamics [1], problems whose mathematical formulation is unknown [2] or real-time simulation performance [3].
However, deep learning is usually very demanding in both data and computational power. A present challenge is not only to develop numerical algorithms that speed-up the already complicated neural network models but also to build smarter architectures which take advantage of the problem structure, potentially reducing the data consumption and improving its generalization. Two historical examples are convolutional [4] and recurrent [5] neural networks, which exploit the structured data in grid elements and timesteps to induce translation and time invariance respectively, both remaining as major breakthroughs in their respective fields.
Motivated by the study of such architectures, which exploit the invariant quantities of the problem, a new machine learning paradigm recently arose known as geometric deep learning [6]. The key insight is to impose specific constraints related to the symmetries of the problem, acting as a strong inductive bias to the learning process. This framework is not only restricted to regular structured data (images or time series) but also to a general case of relations (edges) over arbitrary entities (nodes), mathematically represented with a graph.
In a similar way, physics problems have an intrinsic mathematical structure unveiled by centuries of theoretical and experimental knowledge. From Newton's laws of motion to the Standard Model, the laws of nature can often be described as a set of partial differential equations (PDEs) which can also lead neural networks to find the correct solution of a dynamical problem. These are the foundations of the so-called physics-informed deep learning.

arXiv:2203.01874v3 [cs.LG] 30 May 2022
Thermodynamics-informed Graph Neural Networks The present work aims to combine both inductive biases in order to learn a physical simulator able to predict the dynamics of complex systems in the context of fluid and solid mechanics.

Physics-informed deep learning
Recent works about predicting physics with neural networks [7,1] have demonstrated the convenience of using physical constraints to exploit the dynamical properties of a simulation problem. In fact, the Lagrangian and Hamiltonian reformulations of Newton's classical equations of motion describe a rich mathematical theory based on the Poisson bracket and a symplectic manifold structure. This topological property can be applied to neural networks obtaining robust time integrators, as shown in plenty of recent literature [8,9,10,11,12].
This structure can be extended to non-conservative systems using a dissipative bracket with the so-called GENERIC formalism [13,14]. Several works have already addressed the imposition of this structure bias [15,16,17,3,18], but are hard to learn as the standard multilayer perceptron architectures may not be adequate for larger systems or phenomena whose GENERIC formulation is more complex.

Geometric deep learning
The consideration of a general graph theory framework is motivated by the fact that certain data are not structured in the usual Euclidean space, but in a more complex manifold. For instance, solving PDEs on arbitrary domains is a common problem in physics and engineering, whose solution is rarely found in an Euclidean manifold.
Several graph-based works have achieved great improvements in physics problems such as predicting atomic bond forces [19], particle tracking in high energy physics [20,21], n-body problem with more general interactions [22,23,2], or learning simulators to predict complex fluid interactions [24] and meshed domains [25,26]. These last approaches include the dynamics of the system by predicting the velocity or acceleration, depending on the order of the governing PDE, and integrating them in time. However, they do not include any physical information about the rest of the predicted variables, which remain as black-box direct predictions.
In this work, we propose a learning method of the GENERIC operators and potentials using graph neural networks. This graph-based architecture exploits both the geometrical constraints of discretized systems and the thermodynamically consistent integration scheme of the GENERIC formalism to predict the time evolution of dynamical systems.
The present paper is organized as follows. In section III, we formalize the problem statement and introduce the metriplectic and geometrical inductive biases induced by the proposed architecture. In section IV, we provide three examples of fluid and solid mechanics comparing the results with other methods. Concluding remarks are given in section V.

Problem Statement
Let z ∈ M ⊆ R n be the independent state variables which fully describe a dynamical system up to a certain level of description, with M the set of all the admissible states (space state) which is assumed to have the structure of a differentiable manifold in R n . The given physical phenomenon can be expressed as a system of differential equations encoding the time evolution of its state variables z, where t refers to the time coordinate in the time interval I and F (z, t) refers to an arbitrary nonlinear differential function.
The goal of this paper is to find the convenient mapping F for a dynamical system governed by Eq. (1) from data, in order to efficiently predict its time evolution after a prescribed time horizon T using deep learning [27]. The solution is forced to fulfill the basic thermodynamic requirements of energy conservation and entropy inequality restrictions via the GENERIC formalism.

Metriplectic structure: The GENERIC formalism
In this work, we guarantee the physical meaning of the solution by enforcing the GENERIC structure of the system [13,14], a generalization of the Hamiltonian paradigm to dissipative systems. The GENERIC formulation of time evolution for nonequilibrium systems, described by the set of z state variables, is given by The conservative contribution is assumed to be of Hamiltonian form, requiring an energy potential E = E(z) and a Poisson bracket {z, E} acting on the state vector. This accounts for the reversible dynamics in classical mechanics, expressed by the time structure invariance of the Poisson bracket. Similarly, the remaining irreversible contribution to the energetic balance of the system is generated by the nonequilibrium entropy potential S = S(z) with an irreversible or friction bracket [z, S]. This term is considered an extension of Landau's idea of the time evolution of a state variable towards its equilibrium value via a dissipation potential [28].
For practical use, it is convenient to reformulate the bracket notation using two linear operators: where T * M and T M represent, respectively, the cotangent and tangent bundles of the state space M. These operators inherit the mathematical properties of the original bracket formulation. The operator L = L(z) represents the Poisson bracket and is required to be skew-symmetric (a cosymplectic matrix). Similarly, the friction matrix M = M (z) accounts for the irreversible part of the system and is symmetric and positive semi-definite to ensure the positiveness of the dissipation rate.
Replacing the original bracket formulation in Eq.
(2) with their respective operators, the time-evolution equation for the state variables z is derived, However, the restrictions of L and M are not sufficient to guarantee a thermodynamically consistent description of the dynamics of the system, and two degeneracy conditions are introduced The first one states that the entropy is a degenerate functional of the Poisson bracket, showing the reversible nature of the Hamiltonian contribution to the dynamics. The second expression states the conservation of the total energy of the system with a degenerate condition of the energy with respect to the friction bracket. These restrictions can be reformulated in a matrix form in terms of the L and M operators, resulting in the following degeneracy restrictions: The degeneracy conditions, in addition to the non-negativeness of the irreversible bracket, guarantees the first (energy conservation) and second (entropy inequality) laws of thermodynamics, The specific topological requirements based on a pair of skew-symmetric (L) and symmetric (M ) operators over a smooth manifold is called a metriplectic structure [29,30]. We make use of this hard constraint to construct a thermodynamically-sound integrator, acting as the first inductive bias of our approach. This integration is performed using a forward Euler scheme with time increment ∆t and the GENERIC formalism in Eq. (4), resulting in the following expression In order to learn the GENERIC operators L, M and potentials E, S for each particle of the domain, we propose the use of graph-based deep learning, which exploits the geometrical structure of that specific domain.

Geometric structure: Graph Neural Networks
Let G = (V, E, u) be a directed graph, where V = {1, ..., n} is a set of |V| = n vertices, E ⊆ V × V is a set of |E| = e edges and u ∈ R Fg is the global feature vector. Each vertex and edge in the graph is associated with a node and a pairwise interaction between nodes respectively in a discretized physical system. The global feature vector defines the properties shared by all the nodes in the graph, such as gravity or elastic properties. For each vertex i ∈ V we associate a feature vector v i ∈ R Fv , which represents the physical properties of each individual node. Similarly, for each edge (i, j) ∈ E we associate an edge feature vector e ij ∈ R Fe . In practice, the positional state variables of the system (q i ) are assigned to the edge feature vector e ij so the edge features represent relative distances (q ij = q i − q j ) between nodes, giving a distance-based attentional flavour to the graph network [31,32,33] and translational invariance [34,35]. The rest of the state variables are assigned to the node feature vector v i . The external interactions, such as forces applied to the system, are included in an external load vector f i . A simplified scheme of the graph codification of a physical system is depicted in Fig. 1. Figure 1: Physical system domain discretized in a mesh with node state variables z i , relative nodal distances q ij , external interactions f i and global properties u (top). Graph representation of the same system, with node and edge attributes: v i and e ij (bottom).
These features are fed into an encode-process-decode scheme [23], which consists on several multilayer perceptrons (MLPs) shared between all the nodes and edges of the graph. The algorithm consists of five steps (Fig. 2):

Encoding
We use two MLPs (ε v , ε e ) to transform the vertex and edge initial feature vectors into higher-dimensional embeddings

Processing
The processor is the core task of the algorithm, as it shares the nodal information between vertices via message passing and modifies the hidden vectors in order to extract the desired output of the system. First, a MLP (π e ) computes the updated edge features x ij for each graph edge, based on the current edge features, global features, and sending and recieving node, Then, for each node the messages are pooled with a permutation invariant function φ based on the neighborhood N i = {j ∈ V|(i, j) ∈ E} of the node i. Last, the node embeddings are updated with a second MLP (π v ) using the current node features, the pooled messages, the external load vector and the global features, where (·, ·) denotes vector concatenation and x i and x ij are the updated nodal and edge latent vectors.
The processing step is equivalent to the message passing [36] of 1-step adjacent nodes. In order to get the influence of further graph nodes, the process can be recurrently repeated with both shared or unshared parameters in M processing blocks and optionally using residual connections [37]. In this approach, we use both unshared parameters and residual connections to each message passing block and sum as aggregation function φ. Note that the computed messages φ(x ij ) represent a hidden embedding of the intermolecular interactions of the system (internal messages) whereas the vector f i accounts for the external interactions (external messages).

Decoding
The last block extracts the relevant physical output information y i ∈ R Fy of the system from the node latent feature vector, implemented with a MLP (δ v ). In this work, we predict for each particle the GENERIC energy E and entropy S potentials and the flattened operators l and m:

Reparametrization
A last processing step is needed to get the GENERIC parameters before integrating the state variables. Both operators in matrix form L and M are constructed using the flattened output of the Graph Neural Network l and m respectively, reshaped in lower-triangular matrices. The skew-symmetric and positive semi-definite conditions are imposed by construction using the following parametrization: Both E and S are directly predicted for every node. Then, these potentials can be differentiated with respect to the network input in order to get the gradients ∂E ∂z and ∂S ∂z needed for the GENERIC integrator. These gradients are easily obtained using automatic differentiation [38], and ensures the integrability of the energy and entropy gradients [39].
Considering the dimensions of the lower triangular matrices and the scalar value of both potentials, the output dimension of the decoder network is F y = n(n + 1) 2 + n(n − 1) 2 where n represent the dimension of the state variables z.

Integration
The single-step integration of the state variables of the system z t → z t+1 is then performed using Eq. (8).

Learning procedure
The complete dataset D is composed by N sim multiparametric simulation cases of a dynamical system evolving in time. Each case D i contains the labelled pair of a single-step state vector z t and its evolution in time z t+1 for each node of the system where the dataset D is disjointly partitioned in 80% training, 10% test and 10% validation sets: D train , D val and D test respectively.
The training is performed in a single-snapshot supervision, which has two main advantages: (i) enables parallelization between snapshots, which decreases training time, and (ii) avoids intensive memory usage due to a several-snapshot recursive training. The loss function is divided in two terms: • Data loss: This term accounts for the correct prediction of the state vector time evolution using the GENERIC integrator. It is defined as the MSE along the graph nodes and state variables between the predicted and the ground-truth time derivative of the state vector in a given snapshot, where · 2 denotes the L2-norm. The choice of the time derivative instead of the state vector itself is to regularize the global loss function to a uniform order of magnitude with respect to the degeneracy terms, as shown in Eq. (8).
• Degeneracy loss: This condition is added to the optimization in order to force the degeneracy conditions of the Poisson and dissipative operators, which ensure thermodynamical consistency of the integrator. It is defined as the MSE along the graph nodes and state variables of two residual terms corresponding to the energy and entropy degeneracy conditions, Alternative approaches are found in the literature to impose this degeneracy restrictions, such as a specific tensor parametrization of the brackets [17] or forcing ortogonality using additional skew-symmetric matrices [18]. However, we decide to include it as a soft constraint in order to allow more flexibility in the learning process and improve convergence while maintaining the degeneracy conditions up to an admissible error.
The global loss term is a weighted mean of the two terms over the shuffled N batch batched snapshots, As the energy and entropy are supervised only by their gradients, we remark that (i) they are learnt up to an integration constant value and (ii) the activation functions must have a sufficient degree of continuity. To meet this second requirement, one must select activations with non-zero second derivative in order to have a correct backpropagation of the weights and biases. Thus, linear or rectified units (ReLU, Leaky ReLU, RReLU) are not appropriate for this task. It is well known [40] that logistic functions such as sigmoid and hyperbolic tangent are universal approximators of any derivative arbitrarily well, but are not optimal for very deep neural networks architectures, as they suffer from several problems such as vanishing gradients. Then, the correct activation functions suitable for learning gradients are the ones which combine both non-zero second derivatives and ReLU-type non-linearities, such as Softplus, Swish [41] or Mish [42]. In the present work we use the Swish activation function.
The inputs and outputs of the networks are standardized using the training dataset statistics. Gaussian noise is also added to the inputs during training in order to model the accumulation of error during the time integration [25], which is not contemplated in a single-snapshot training, with the variance of the noise σ 2 noise as a tunable hyperparameter and zero mean value. All the cases are optimized using Adam [43] and a multistep learning rate scheduler.
The code is fully implemented in Pytorch. Our datasets and trained networks are publicly available online at https: //github.com/quercushernandez.

Evaluation metrics
Two ablation studies are performed to evaluate the method presented in this work. The first case is performed using only Graph Neural Networks (from now, GNN) with similar architecture and learning procedure used in prior works [24,25] and no metriplectic integrator. In the second case, we impose the metriplectic structure [15,16] (from now, SPNN), using standard MLPs with no graph computations. Both alternative methods are tuned for equal parameter count in order to get a fair comparison of the results.
All the results are computed with the integration scheme in Eq. (8) iteratively from the initial conditions to the prescribed time horizon T , denoted as rollout. The rollout prediction error is quantified by the relative L2 error, computed with Eq. (19) for each snapshot and simulation case, The results are represented in Fig. 9, 10 and 11 showing the rollout statistics for all the snapshots divided in train and test simulations, state variables and method used (Ours, GNN or SPNN).

Description
The first example is a shear (Couette) flow of an Oldroyd-B fluid model (Fig. 3). This is a constitutive model for viscoelastic fluids, considering linear elastic dumbbells as a proxy representation of polymeric chains immersed in a solvent.
The edge feature vector contains the relative position of the nodes whereas the rest of the state variables are part of the node feature vector. An additional one-hot vector n is added to the node features in order to represent the boundary and fluid nodes. The global feature vector u represent the Weissenberg and Reynolds numbers of each simulation, resulting in the following feature vectors:

Database and Hyperparameters
The training database for the Couette flow is generated with the CONNFFESSIT technique [44], based on the Fokker-Plank equation [45], using a Monte Carlo algorithm. The fluid is discretized in the vertical direction with N e = 100 elements and N = 101 nodes in a total height of H = 1. A total of 10,000 dumbells are considered at each nodal location in the model. The lid velocity is set to V = 1, with variable Weissenberg We ∈ [1,2] and Reynolds number Re ∈ [0.1, 1], summing a total of N sim = 100 cases. The simulation is discretized in N T = 150 time increments of ∆t = 6.7 · 10 −3 . Following Eq. (21), the dimensions of the graph feature vectors are F e = 3, F v = 5 and F g = 2. The hidden dimension of the node and edge latent vectors is F h = 10. The learning rate is set to l r = 10 −3 with decreasing order of magnitude on epochs 2000 and 4000, and a total number of N epoch = 6000. The training noise variance is set to σ 2 noise = 10 −2 .

Results
The rollout results for the Couette flow are presented in Fig. 9. A substantial improvement is shown in the present approach over the two other methods, which remain in a similar performance. Note that the skewed distributions towards higher errors on each box is due to the error accumulation on snapshots further in time from the starting conditions, where errors are lower. Fig. 8 (left) shows that the degeneracy conditions imposed by our method ensure the thermodynamical consistency of the learnt energy and entropy potentials.

Description
The next example is a viscoelastic cantilever beam subjected by a bending force. The material is characterized by a single-term polynomial strain energy potential, described by the following equation where U is the strain energy potential, J el is the elastic volume ratio, I 1 and I 2 are the two invariants of the left Cauchy-Green deformation tensor, C 10 and C 01 are shear material constants and D 1 is the material compressibility parameter. The viscoelastic component is described by a two-term Prony series of the dimensionless shear relaxation modulus, with relaxation coefficients ofḡ 1 andḡ 2 , and relaxation times of τ 1 and τ 2 . The state variables for the viscoelastic beam on each node are the position q, velocity v and stress tensor σ, The relative deformed position is included into the edge feature vector whereas the rest of the variables are part of the node feature vector. An additional one-hot vector n is added to the node features in order to represent the encastre and beam nodes. The external load vector F is included in the node processor MLP as an external interaction. No global feature vector is needed in this case, resulting in the following feature vectors:
Following Eq. (25), the dimensions of the graph feature vectors are F e = 4, F v = 11 and F g = 0. The hidden dimension of the node and edge latent vectors is F h = 50. The learning rate is set to l r = 10 −4 with decreasing order of magnitude on epochs 600 and 1200, and a total number of N epoch = 1800. The training noise variance is set to σ 2 noise = 10 −5 .

Results
The rollout results for the bending viscoelastic beam are presented in Fig. 10. The errors achieved by the present approach are again below the other two methods. The beam deformed configuration of three different test simulation snapshots are represented in Fig. 5, with the color code representing the xx component of the stress tensor. Similarly to the previous case, Fig. 8 (center) shows the thermodynamical consistency of our dynamical integration. The state variables for the flow past a cylinder are the velocity v and the pressure field P , The flow is computed with an Eulerian description of the output fields. Thus, the nodal coordinates (q 0 ) are fixed in space and considered as edge features, whereas the whole state variables are assigned to the node features. An additional one-hot vector n is added to the node features in order to represent the inlet/outlet, walls or fluid nodes. No global feature vector is needed in this case, resulting in the following feature vectors:

Database and Hyperparameters
The

Results
The rollout results for the flow past a cylinder problem are presented in Fig. 11. In this example the domain varies significantly, using a different unstructured mesh for each simulation. Thus, the graph-based architectures outperform the vanilla SPNN, which is meant for fixed structured problems. Considering the other two methods, our approach outperforms the standard GNN architecture due to the metriplectic structure imposition over the dynamical problem, as depicted in Fig. 8. A single snapshot of the whole rollout of three different test simulations are represented in Fig. 7, with the color code representing the x component of the velocity field.

Conclusions
We have presented a method to predict the time evolution of an arbitrary dynamical system based on two inductive biases. The metriplectic bias ensures the correct thermodynamic structure of the integrator based on the GENERIC formalism, whose operators and potentials are estimated using computations over graphs, i.e. exploiting the geometric structure of the problem. The results show relative mean errors of less than 3% in all the tested examples, outperforming two other state-of-the-art techniques based on only physics-informed and geometric deep learning respectively. These results confirm that both biases are necessary to achieve higher precision in the predicted simulations. The use of both techniques combine the computational power of geometric deep learning with the rigorous foundation of the GENERIC formalism, which ensure the thermodynamical consistency of the results.
The limitations of the presented technique are related to the computational complexity of the model. Large simulations with fine grids require a high amount of message passing to get the information across the whole domain, or a very Future work may overcome the stated limitations by combining graph representations with model order reduction techniques, such as autoencoders or U-net architectures [46,47]. The idea is to replace deep message passing with various coarse-graining steps, allowing the boundary information to reach every node in the simulation domain while reducing the number of parameters of the neural network. Another interesting topic to extend our work is to improve the generalization and decrease the amount of training data via equivariant arquitectures [48,49], which avoid data augmentation by exploiting the invariance to certain groups such as rotations SO(3) or general Euclidean transformations E(3). As the present work is only limited to in-silico experiments, future work may extend the proposed method to measured datasets in real-world applications, such as digital twins of industrial processes or real-time augmented/virtual reality environments.