Quercus Hernández
email: quercus@unizar.es

Alberto Badías
email: abadias@unizar.es

David González

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Elías Cueto
email: ecueto@unizar.es

Structure-preserving neural networks ✩

Keywords:

We develop a method to learn physical systems from data that employs feedforward neural networks and whose predictions comply with the first and second principles of thermodynamics. The method employs a minimum amount of data by enforcing the metriplectic structure of dissipative Hamiltonian systems in the form of the socalled General Equation for the Non-Equilibrium Reversible-Irreversible Coupling, GENERIC (Öttinger and Grmela (1997) [36]). The method does not need to enforce any kind of balance equation, and thus no previous knowledge on the nature of the system is needed. Conservation of energy and dissipation of entropy in the prediction of previously unseen situations arise as a natural by-product of the structure of the method. Examples of the performance of the method are shown that comprise conservative as well as dissipative systems, discrete as well as continuous ones.

Introduction

With the irruption of the so-called fourth paradigm of science [START_REF] Hey | The Fourth Paradigm: Data-Intensive Scientific Discovery[END_REF] a growing interest is detected on the machine learning of scientific laws. A plethora of methods have been developed that are able to produce more or less accurate predictions about the response of physical systems in previously unseen situations by employing techniques ranging from classical regression to the most sophisticated deep learning methods.

For instance, recent works in solid mechanics have substituted the constitutive equations with experimental data [START_REF] Ayensa-Jiménez | A new reliability-based data-driven approach for noisy experimental data with physical constraints[END_REF][START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF], while conserving the traditional approach on physical laws with high epistemic value (i.e., balance equations, equilibrium). Similar approaches have applied this concept to the unveiling (or correction) of plasticity models [START_REF] Ibáñez | Hybrid constitutive modeling: datadriven learning of corrections to plasticity models[END_REF], while others created the new concept of constitutive manifold [START_REF] Ibanez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF][START_REF] Ibañez | Datadriven non-linear elasticity: constitutive manifold construction and problem discretization[END_REF]. Other approaches are designed to unveil an explicit, closed form expression for the physical law governing the phenomenon at hand [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF].

An interest is observed in the incorporation of the already existing scientific knowledge to these data-driven procedures. This interest is two-fold. Indeed, we prefer not to get rid of centuries of scientific knowledge and rely exclusively on powerful machine learning strategies. Existing theories have proved to be useful in the prediction of physical phenomena and are still in the position of helping to produce very accurate predictions. This is the procedure followed in the so-called data-driven computational mechanics approach mentioned before. On the other hand, these theories help to keep the consumption of data to a minimum. Data are expensive to produce and to maintain. Already existing scientific knowledge could alleviate the amount of data needed to produce a successful prediction.

The mentioned works on data-driven computational mechanics usually rely on traditional machine learning algorithms, which are very precise and tested but usually computationally expensive. With the recent advances in data processing, computing resources and machine learning, neural networks have become a powerful tool to analyze traditionally hard problems such as image classification [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Shin | Deep convolutional neural networks for computer-aided detection: Cnn architectures, dataset characteristics and transfer learning[END_REF], speech recognition [START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF][START_REF] Hinton | Deep neural networks for acoustic modeling in speech recognition: the shared views of four research groups[END_REF] or data compressing [START_REF] Romero | Quantum autoencoders for efficient compression of quantum data[END_REF][START_REF] Theis | Lossy image compression with compressive autoencoders[END_REF]. These new machine learning methods outperform many of the traditional ones, both in modeling capacity and computational time (once trained, certain neural networks can easily handle real time requirements). Recent works in the machine learning community [START_REF] Yoon | Gradient-based inference for networks with output constraints[END_REF][START_REF] Márquez-Neila | Imposing hard constraints on deep networks: promises and limitations[END_REF][START_REF] Nandwani | A primal dual formulation for deep learning with constraints[END_REF] have shown that neural networks are also versatile in constraint optimizations. This is the approach followed by several authors in the context of physical simulations, which aim to solve a set of partial differential equations (PDEs) in complex dynamical systems. Physical problems must satisfy inherently certain conditions dictated by physics, often formulated as conservation laws, and can be imposed to a neural network using extra loss terms in the constrained optimization process [START_REF] Magiera | Constraint-aware neural networks for Riemann problems[END_REF].

Similar constraints are imposed in the so-called physically-informed neural networks approach [START_REF] Raissi | Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Zhang | Learning in modal space: solving time-dependent stochastic pdes using physics-informed neural networks[END_REF]. This family of methods employs neural networks to solve highly nonlinear partial differential equations (PDEs) resulting in very accurate and numerically stable results. However, they rely on prior knowledge of the governing equations of the problem.

The authors have introduced the so-called thermodynamically consistent data-driven computational mechanics [START_REF] Ghnatios | Data-driven generic modeling of poroviscoelastic materials[END_REF][START_REF] González | Consistent data-driven computational mechanics[END_REF][START_REF] González | Learning corrections for hyperelastic models from data[END_REF]. Unlike other existing works, this approach does not impose any particular balance equation to solve for. Instead, it relies on the imposition of the right thermodynamic structure of the resulting predictions, as dictated by the so-called GENERIC formalism [START_REF] Grmela | Dynamics and thermodynamics of complex fluids. I. Development of a general formalism[END_REF]. As will be seen, this ensures conservation of energy and the right amount of entropy dissipation, thus giving rise to predictions satisfying the first and second principles of thermodynamics. These techniques, however, employ regression to unveil the thermodynamic structure of the problem at the sampling points. For previously unseen situations, they employ interpolation on the matrix manifold describing the system.

Recent works in symplectic networks [START_REF] Jin | Symplectic networks: intrinsic structure-preserving networks for identifying Hamiltonian systems[END_REF] have by-passed those drawbacks by exploiting the mathematical properties of Hamiltonian systems, so no prior knowledge of the system is required. However, this technique only operates on conservative systems with no entropy generation.

The aim of this work is the development of a new structure-preserving neural network architecture capable of predicting the time evolution of a system based on experimental observations on the system, with no prior knowledge of its governing equations, to be valid for both conservative and dissipative systems. The key idea is to merge the proven computational power of neural networks in highly nonlinear physics with thermodynamic consistent data-driven algorithms. The resulting methodology, as will be seen, is a powerful neural network architecture, conceptually very simple-based on standard feedforward methodologies-that exploits the right thermodynamic structure of the system as unveiled from experimental data, and that produces interpretable results [START_REF] Murdoch | Interpretable machine learning: definitions, methods, and applications[END_REF].

The outline of the paper is as follows. A brief description of the problem setup is presented in Section 2. Next, in Section 3, the methodology is presented of both the GENERIC formalism and the feed-forward neural networks used to solve the stated problem. This technique is used in different physical systems of increasing complexity: a double thermoelastic pendulum (Section 4) and a Couette flow in a viscoelastic fluid (Section 5). The paper is completed with a discussion in Section 6.

Problem statement

Weinan E seems to be the first author in interpreting the process of learning physical systems as the solution of a dynamical system [START_REF] Weinan | A proposal on machine learning via dynamical systems[END_REF]. Consider a system whose governing variables will be hereafter denoted by z ∈ M ⊆ R n , with M the state space of these variables, which is assumed to have the structure of a differentiable manifold in R n .

The problem of learning a given physical phenomenon can thus be seen as the one of finding an expression for the time evolution of their governing variables z,

ż = dz dt = F (x, z, t), x ∈ ∈ R D , t ∈ I = (0, T ], z(0) = z 0 , (1)
where x and t refer to the space and time coordinates within a domain with D = 2, 3 dimensions. F (x, z, t) is the function that gives, after a prescribed time horizon T , the flow map z 0 → z(z 0 , T ).

While this problem can be seen as a general supervised learning problem (we fix both z 0 and z), when we have additional information about the physics being represented by the sought function F , it is legitimate to try to include it in the search procedure. W. E seems to have been the first in suggesting to impose a Hamiltonian structure on F if we know that energy is conserved, for instance [START_REF] Weinan | A proposal on machine learning via dynamical systems[END_REF]. Very recently, two different approaches follow this same rationale [START_REF] Bertalan | On learning Hamiltonian systems from data[END_REF][START_REF] Jin | Symplectic networks: intrinsic structure-preserving networks for identifying Hamiltonian systems[END_REF].

For conservative systems, therefore, imposing a Hamiltonian structure seems a very appealing way to obtain thermodynamics-aware results. However, when the system is dissipative, this method does not provide with valid results. Given the importante of dissipative phenomena (viscous solids, fluid dynamics, ...) we explore the right thermodynamic structure to impose to the search methodology.

The goal of this paper is to develop a new method of solving Eq. (1) using state of the art deep learning tools, in order to predict the time evolution of the state variables of a given system. The solution is forced to fulfill the basic thermodynamic requirements of energy conservation and entropy inequality restrictions via the GENERIC formalism, presented in the next section.

Methodology

In this section we develop the appropriate thermodynamic structure for dissipative systems. Classical systems modeling can be done at a variety of scales. We could think of the most detailed (yet often impractical) scale of molecular dynamics, where energy conservation applies and the Hamiltonian paradigm can be imposed. However, the number of degrees of freedom and, noteworthy, the time scale, renders this approach of little interest for many applications. On the other side of the spectrum lies thermodynamics, where only conserved, invariant, quantities are described and thus there is no need for conservation principles. At any other (mesoscopic) scale, unresolved degrees of freedom give rise to the appearance of fluctuation in the results (or its equivalent, dissipation). At these scales, traditional modeling procedures imply expressing physical insights in the form of governing equations [START_REF] Grmela | Generic guide to the multiscale dynamics and thermodynamics[END_REF]. These equations are then validated from experimental observations.

Alternatively, thermodynamics can be thought of as a meta-physics, in the sense that it is actually a theory of theories [START_REF] Grmela | Gradient and generic evolution towards reduced dynamics[END_REF]. It provides us with the right theoretic framework in which basic principles are met. And, in particular for any of these intermediate or mesoscopic scales, a so-called metriplectic structure emerges. The term metriplectic comes for the combination of symplectic and Riemannian (metric) geometry and emphasizes the fact that there are conservative as well as dissipative contributions to the general evolution of such a system. Once such a geometric structure is found for the system, we are in the position of fixing the framework in which our neural networks can look for the adequate prediction of the future states of the system. The particular metriplectic structure that we employ for such a task is known, as stated before, as GENERIC.

The GENERIC formalism

The "General Equation for Non-Equilibrium Reversible-Irreversible Coupling", GENERIC, formalism [START_REF] Grmela | Dynamics and thermodynamics of complex fluids. I. Development of a general formalism[END_REF][START_REF] Hans Christian Öttinger | Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism[END_REF] establishes a mathematical framework in order to model the dynamics of a system. Furthermore, it is compatible with classical equilibrium thermodynamics [START_REF] Christian | Beyond Equilibrium Thermodynamics[END_REF], preserving the symmetries of the system as stated in Noether's theorem. It has served as the basis for the development of several consistent numerical integration algorithms that exploit these desirable properties [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF][START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF].

The GENERIC structure for the evolution in Eq. ( 1) is obtained after finding two algebraic or differential operators

L : T * M → T M, M : T * M → T M,
where T * M and T M represent, respectively, the cotangent and tangent bundles of M. As in general Hamiltonian systems, there will be an energy potential, which we will denote hereafter by E(z). In order to take into account the dissipative effects, a second potential (the so-called Massieu potential) is introduced in the formulation. It is, of course, the entropy potential of the GENERIC formulation, S(z). With all these ingredients, we arrive at a description of the dynamics of the system of the type

dz dt = L ∂ E ∂ z + M ∂ S ∂ z .
(

) 2
As shown in Eq. ( 2), the time evolution of the system described by the nonlinear operator F (x, z, t) presented in Eq. ( 1) is now split in two separated terms:

• Reversible Term: It accounts for all the reversible (non-dissipative) phenomena of the system. In the context of classical mechanics, this term is equivalent to Hamilton's equations of motion that relates the particle position and momentum. The operator L(z) is the Poisson matrix-it defines a Poisson bracket-and is required to be skew-symmetric (a cosymplectic matrix).

• Non-Reversible Term: The rest of the non-reversible (dissipative) phenomena of the system are modeled here. The operator M(z) is the friction matrix and is required to be symmetric and positive semi-definite.

The GENERIC formulation of the problem is completed with the following so-called degeneracy conditions

L ∂ S ∂ z = M ∂ E ∂ z = 0. (3)
The first condition expresses the reversible nature of the L contribution to the dynamics whereas the second requirement expresses the conservation of the total energy by the M contribution. This means no other thing that the energy potential does not contribute to the production of entropy and, conversely, that the entropy functional does not contribute to reversible dynamics. This mutual degeneracy requirement in addition to the already mentioned L and M matrix requirements ensure that:

∂ E ∂t = ∂ E ∂ z • ∂ z ∂t = ∂ E ∂ z L ∂ E ∂ z + M ∂ S ∂ z = 0,
which expresses the conservation of energy in an isolated system, also known as the first law of thermodynamics. Applying the same reasoning to the entropy S:

∂ S ∂t = ∂ S ∂ z • ∂ z ∂t = ∂ S ∂ z L ∂ E ∂ z + M ∂ S ∂ z = ∂ S ∂ z M ∂ S ∂ z ≥ 0,
which guarantees the entropy inequality, this is, the second law of thermodynamics.

Proposed integration algorithm

Once the learning procedure is accomplished, our neural network is expected to integrate the system dynamics in time, given previously unseen initial conditions. In order to numerically solve the GENERIC equation, we formulate the discretized version of Eq. ( 2) following previous works [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF]:

z n+1 -z n t = L • DE Dz + M • DS Dz . ( 4
)
The time derivative of the original equation is discretized with a forward Euler scheme in time increments t, where z n+1 = z t+ t . L and M are the discretized versions of the Poisson and friction matrices. Last, DE Dz and DS Dz represent the discrete gradients, which can be approximated in a finite element sense as:

DE Dz Az, DS Dz B z, (5)
where A and B represent the discrete matrix form of the gradient operators. Finally, manipulating algebraically Eq. ( 4) with Eq. ( 5) and including the degeneracy conditions of Eq. ( 3), the proposed integration scheme for predicting the dynamics of a physical system is the following

z n+1 = z n + t (L • Az n + M • B z n ) (6) subject to: L • B z n = 0, M • Az n = 0,
ensuring the thermodynamical consistency of the resulting model.

To sum up, the main objective of this work is to compute the form of the A(z) and B(z) gradient operator matrices, subject to the degeneracy conditions, in order to integrate the initial system state variables z 0 over certain time steps t of the time interval I. Usually, the form of matrices L and M is known in advance, given the vast literature in the field. If necessary, these terms can also be computed [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF].

Feed-forward neural networks

In the introduction we already mentioned the intrinsic power of neural networks in many fields. The main reason under the fact that neural networks are able to learn and reproduce such a variety of problems is that they are considered to be universal approximators [START_REF] Choi | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], meaning that they are capable of approximating any measurable function to any desired degree of accuracy. The main limitation of this technique is the correct selection of the tuning parameters of the network, also called hyperparameters.

Another universal approximator are polynomials, as they can approximate any infinitely differentiable function as a Taylor power series expansion. The main difference is that neural networks rely on composition of functions rather than sum of power series: [2] • f [1] )(x). [START_REF] Ghnatios | Data-driven generic modeling of poroviscoelastic materials[END_REF] Eq. [START_REF] Ghnatios | Data-driven generic modeling of poroviscoelastic materials[END_REF] shows that the desired output ŷ from a defined input x of a neural network is a composition of different functions f [l] as building blocks of the network in L total layers. The challenge is to select the best combination of functions in the correct order such that it approximates the solution of the studied problem. The simplest building block of artificial deep neural network architectures is the neuron or perceptron (Fig. 1,left). Several neurons are stacked in a multilayer perceptron (MLP), which is mathematically defined as follows x [l] = σ (w [l] x [l-1] + b [l] ), [START_REF] Glorot | Deep sparse rectifier neural networks[END_REF] where l is the index of the current layer, x [l-1] and x [l] are the layer input and output vector respectively, w [l] is the weight matrix of the last layer, b [l] is the bias vector of the last layer and σ is the activation function. If no activation function is applied, the MLP is equivalent to a linear operator. However, σ is chosen to be a nonlinear function in order to increase the capacity of modeling more complex problems, which are commonly nonlinear. In classification problems, the traditional activation function is the logistic function (sigmoid) whereas in regression problems, Rectified Linear Unit (ReLU) [START_REF] Glorot | Deep sparse rectifier neural networks[END_REF] or hyperbolic tangent are commonly used.

ŷ = ( f [L] • f [L-1] • ... • f [l] • ... • f
In this work, we use a deep neural network architecture known as feed-forward neural network [START_REF] Schmidhuber | Deep learning in neural networks: an overview[END_REF]. It consists of a several layer of multilayer perceptrons with no cyclic connections, as shown in Fig. 1

(right).

The input of the neural net is the vector state of a given time step z n , and the outputs are the concatenated GENERIC matrices A net n and B net n : for a system with n state variables the number of inputs and outputs are N in = n and N out = 2n 2 .

Then, using the GENERIC integration scheme, the state vector at the next time step z net n+1 is obtained. This method is repeated for the whole simulation time T with a total of N T snapshots.

The state variables of a general dynamical system may differ in several orders of magnitude from each other, due to their own physical nature or measurement units. Then, a pre-processing of the input data (scaling or normalization) can improve the model performance and stability.

The number of hidden layers N h depends on the complexity of the problem. Increasing the net size raises the computational power of the net to model more complex phenomena. However, it slows the training process and could lead to data overfitting, limiting its generalization and extrapolation capacity. The size of the hidden layers is chosen to be the same as the output size of the net N out .

The cost function for our neural network is composed of three different terms:

• Data loss: The main loss condition is the agreement between the network output and the real data. It is computed as the squared error sum, computed between the predicted state vector z net n+1 and the ground truth solution z GT n+1 for each time step.

L data n = z GT n+1 -z net n+1 2 2 . ( 9
)
• Fulfillment of the degeneracy conditions: The cost function will also account for the degeneracy conditions in order to ensure the thermodynamic consistency of the solution, implemented as the sum of the squared elements of the degeneracy vectors for each time step,

L degen n = L • B net n z net n 2 2 + M • A net n z net n 2 2 . ( 10
)
This term acts as a regularization of the loss function and, at the same time, is the responsible of ensuring thermodynamic consistency. So to speak, it is the cornerstone of our method.

• Regularization: In order to avoid overfitting, an extra L2 regularization term L reg is added to the loss function, defined as the sum over the squared weight parameters of the network.

L reg = L l n [l] i n [l+1] j (w [l] i, j ) 2 . ( 11
)
The total cost function is computed as the sum squared error (SSE) of the data loss and degeneracy residual, in addition to the regularization term, at the end of the simulation time T for each train case. The regularization loss is highly dependent on the size of the network layers and has different scaling with respect to the other terms, so it is compensated with the regularization hyperparameter (weight decay) λ r . An additional weight λ d is added to the data loss term, which accounts for the relative scaling error with respect to the degeneracy conditions.

L = N T n=0 (λ d L data n + L degen n ) + λ r L reg . ( 12
)
The usual backpropagation algorithm [START_REF] Paszke | Automatic differentiation in pytorch[END_REF] is then used to calculate the gradient of the loss function for each net parameter (weight and bias vectors), which are updated with the gradient descent technique [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF]. The process is then repeated for a maximum number of epochs n epoch . The resulting training algorithm is sketched in Fig. 2.

The proposed methodology is tested with two different databases of nonlinear physical systems, split in a partition of train cases (N train = 80% of the database) and test cases (N test = 20% of the database). The net performance is evaluated with the mean squared error (MSE) of the state variables prediction, associated with the data loss term, Eq. ( 9), over all the time snapshots,

MSE data (z i ) = 1 N T N T n=0 z GT i,n -z net i,n 2 . ( 13
)
The same procedure is applied to the degeneracy constraint, associated with the degeneracy loss term, Eq. [START_REF] González | Learning corrections for hyperelastic models from data[END_REF], over all the time snapshots,

MSE degen (z i ) = 1 N T N T n=0 L • B net i,n z net i,n + M • A net i,n z net i,n . ( 14
)
As a general error magnitude of the algorithm, the average MSE of both the train (N = N train ) and test trajectories (N = N test ) is also reported for both the data (m = data) and degeneracy (m = degen) constraints, MSE

m (z) = 1 N N i=1 MSE m (z i ). ( 15
)
Algorithm 1 and Algorithm 2 show a pseudocode of our proposed algorithm to both the training and test processes. The proposed method is fully implemented in PyTorch [START_REF] Paszke | Pytorch: an imperative style, high-performance deep learning library[END_REF] and trained in an Intel Core i7-8665U CPU.

Validation examples: double thermo-elastic pendulum

Description

The first example is a double thermo-elastic pendulum (Fig. 3) consisting of two masses m 1 and m 2 connected by two springs of variable lengths λ 1 and λ 2 and natural lengths at rest λ 0 1 and λ 0 2 . The set of variables describing the double pendulum are here chosen to be

S = {z = (q 1 , q 2 , p 1 , p 2 , s 1 , s 2 ) ∈ (R 2 × R 2 × R 2 × R 2 × R × R), q 1 = 0, q 1 = q 2 }, (16)
where q i , p i and s i are the position, linear momentum and entropy of each mass i = 1, 2. Initialize losses: L data , L degen = 0;

for snapshot ← 1, N T do Forward propagation: [A net n , B net n ] ← Net(z GT n );
Eq. ( 8) Time integration:

z net n+1 ← z net n + t (L • A net n z net n + M • B net n z net n );
Eq. ( 4)

Update data loss:

L data ← L data + L data n ;
Eq. ( 9)

Update degeneracy loss:

L degen ← L degen + L degen n ;
Eq. ( 10) end for

SSE loss function: L ← λ d L data + L degen + λ r L reg
Eq. ( 11), Eq. ( 12

z net 0 ← z GT 0 ; for snapshot ← 1, N T do Forward propagation: [A net n , B net n ] ← Net(z net n );
Eq. ( 8) Time step integration:

z net n+1 ← z net n + t (L • A net n z net n + M • B net n z net n );
Eq. (4) Update state vector:

z net n ← z net n+1 ; Update snapshot: n ← n + 1;
end for Compute MSE data , MSE degen ; Eq. ( 13), Eq. ( 14) end for Compute MSE data , MSE degen ; Eq. ( 15) Fig. 3. Double thermo-elastic pendulum.

The lengths of the springs λ 1 and λ 2 are defined solely in terms of the positions as

λ 1 = √ q 1 • q 1 , λ 2 = (q 2 -q 1 ) • (q 2 -q 1 ).
The total energy of the system can be expressed as the sum of the kinetic energy of the two masses K i and the internal energy of the springs e i for i = 1, 2, The total entropy of the double pendulum is the sum of the entropies of the two masses s i ,

E = E(z) = i K i (z) + i e i (λ i , s i ), K i = 1 2m i |p i | 2 . ( 17
)
S = S(z) = s 1 + s 2 . ( 18
)
This model includes thermal effects in the stretching of the springs due to the Gough-Joule effect. The absolute temperatures T i at each spring are obtained through Eq. ( 19). These temperature changes induce a heat flux between both springs, being proportional to the temperature difference and a conductivity constant κ > 0,

T i = ∂e i ∂s i . ( 19
)
In this case, there is a clear contribution of both conservative Hamiltonian mechanics (mass movement) and non-Hamiltonian dissipative effects (heat flux), resulting in a non-zero Poisson matrix (M = 0). Thus, the GENERIC matrices associated with this physical system are known to be [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF] L = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -1/2 0 0 0 0 -1/2 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 20
)

Database and hyperparameters

The training database is generated with a thermodynamically consistent time-stepping algorithm [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF] in MATLAB. The masses of the double pendulum are set to m 1 = 1 kg and m 2 = 2 kg, joint with springs of a natural length of λ 0 1 = 2 m and λ 0 2 = 1 m and thermal constant of C 1 = 0.02 J and C 2 = 0.2 J and conductivity constant of κ = 0.5. The simulation time of the movement is T = 60 s in time increments of t = 0.3 s (N T = 200 snapshots).

The database consists of the state vector, Eq. ( 16), of 50 different trajectories with random initial conditions of position q i and linear momentum p i of both masses m i (i = 1, 2) around a mean position and linear momentum of q 1 = [4.5, 4.5] m, p 1 = [2, 4.5] kg•m/s, and q 2 = [-0.5, 1.5] m, p 2 = [1.4, -0.2] kg•m/s respectively. Although the initial conditions of the simulations are similar, it results in a wide variety of the mass trajectories due to the chaotic behavior of the system. This database is split randomly in 40 train trajectories and 10 test trajectories. Thus, there is a total of 80.000 training snapshots and 20.000 test snapshots.

The net input and output size is N in = 10 and N out = 2N 2 in = 200. The state vector is normalized based on the training set statistical mean and standard deviation. The number of hidden layers is N h = 5 with ReLU activation functions and linear in the last layer. It is initialized according to the Kaiming method [START_REF] He | Delving deep into rectifiers: surpassing human-level performance on imagenet classification[END_REF] with normal distribution and the optimizer used is Adam [START_REF] Diederik | Adam: a method for stochastic optimization[END_REF] with a weight decay of λ r = 10 . Time evolution of the state variables in a test trajectory of a double thermo-elastic pendulum using a time-stepping solver (Ground Truth, GT) and the proposed GENERIC integration scheme (SPNN). Since every variable has a vectorial character, both components are depicted and labeled as X and Y , respectively.

Results

Fig. 5 shows the time evolution of the state variables (position, momentum and entropy) of each mass given by the solver and the neural net.

Table 1 shows the mean squared error of the data and degeneration loss terms for all the state variables of the double pendulum. The results are computed separately as the mean over all the train and test trajectories using Eq. ( 15). Fig. 6 and Fig. 7 show the time evolution of the internal and kinetic energy (Eq. ( 17)) and the entropy (Eq. ( 18)) respectively for the two pendulum masses (i = 1, 2). The total energy is conserved and the total entropy satisfies the entropy inequality, fulfilling the first and second laws of thermodynamics respectively. The mean error for both train and test trajectories is reported in Table 2.

Fig. 6. Time evolution of the energy in a test trajectory of a double thermo-elastic pendulum using a time-stepping solver (Ground Truth, GT) and the proposed GENERIC integration scheme (Net). Fig. 7. Time evolution of the entropy in a test trajectory of a double thermo-elastic pendulum using a time-stepping solver (Ground Truth, GT) and the proposed GENERIC integration scheme (SPNN).

Couette flow of an Oldroyd-B fluid

Description

The second example is a shear (Couette) flow of an Oldroyd-B fluid model. This is a constitutive model for viscoelastic fluids, consisting of linear elastic dumbbells (representing polymer chains) immersed in a solvent.

The Oldroyd-B model arises in the modeling of flows of diluted polymeric solutions. This model can be obtained both from a purely macroscopic point of view as well as from a microscopic one, by modeling polymer chains as linear dumbbells diluted in a Newtonian substrate. Alternatively, it can also be obtained by considering the deviatoric part T of the stress tensor σ (the so-called extra-stress tensor), to be of the form where the triangle denotes the non-linear Oldroyd's upper-convected derivative [START_REF] Owens | Computational Rheology[END_REF]. Coefficients η 0 , λ 1 and λ 2 are model parameters. It is standard to denote the strain rate tensor by γ = (∇ s v) = D.

T + λ 1 ∇ T = η 0 γ + λ 2 ∇ γ , (21)
Finally, the stress in the solvent (denoted by a subscript s) and polymer (denoted by a subscript p) are given by

T = η s γ + τ , so that τ + λ 1 ∇ τ = η p γ ,
which is the constitutive equation for the elastic stress.

Pseudo-experimental data are obtained by the CONNFFESSIT technique [START_REF] Laso | Calculation of viscoelastic flow using molecular models: the connffessit approach[END_REF], based on the Fokker-Plank equation [START_REF] Le Bris | Multiscale modelling of complex fluids: a mathematical initiation[END_REF]. This equation is solved by converting it in its corresponding Itô stochastic differential equation,

dr x = ∂ v ∂ y r y - 1 2We r x dt + 1 √ We dV t , dr y = - 1 2We r y dt + 1 √ We dW t , ( 22
)
where v is the flow velocity, r = [r x , r y ] , r x = r x (y, t) the position vector and assuming a Couette flow so that r y = r y (t)

depends only on time, we stand for the Weissenberg number and V t , W t are two independent one-dimensional Brownian motions. This equation is solved via Monte Carlo techniques, by replacing the mathematical expectation by the empirical mean.

The model relies on the microscopic description of the state of the dumbbells. Thus, it is particularly useful to base the microscopic description on the evolution of the conformation tensor c = rr , this is, the second moment of the dumbbell end-to-end distance distribution function. This tensor is in general not experimentally measurable and plays the role of an internal variable. The expected xy stress component tensor will be given by

τ = We 1 K K k=1 r x r y ,
where K is the number of simulated dumbbells and =

ν p ν p
is the ratio of the polymer to solvent viscosities.

The state variables selected for this problem are the position of the fluid on each node of the mesh, see Fig. 8, its velocity v in the x direction, internal energy e and the conformation tensor shear component τ ,

S = {z = (q, v, e, τ ) ∈ (R 2 × R × R × R)}. (23)
The GENERIC matrices associated with each node of this physical system are the following

L = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 0 0 0 0 0 0 -1 0 0 0 1 -1 0 0 -1 0 0 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ , M = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎦ . ( 24
)
In order to simulate a measurement of real captured data, Gaussian noise is added to the state vector, computed as a random variable following a normal distribution with zero mean and standard deviation proportional to the standard deviation of the database σ z and noise level ν,

z GT noise = z GT + ν • σ z • N (0, 1) (25)
The results of both the noise-free and the noisy database are compared with two different network architectures:

• Unconstrained network: This architecture is the same as the proposed network but removing the degeneracy conditions of the energy and entropy, Eq. ( 10), in the loss function. These conditions ensure the thermodynamic consistency of the resulting integrator, so not including them affects negatively in the accuracy of the results, as will be seen.

• Black-Box network: In this case, no GENERIC architecture is imposed, acting as a black-box integrator trained to directly predict the state vector time evolution z t+1 from the previous time step z t . This naive approach is shown to be inappropriate, as no physical restrictions are given to the model.

Database and hyperparameters

The training database for this Oldroyd-B model is generated in MATLAB with a multiscale approach [START_REF] Le Bris | Multiscale modelling of complex fluids: a mathematical initiation[END_REF] in the dimensionless form. The fluid is discretized in the vertical direction with N = 100 elements (101 nodes) in a total height of H = 1. A total of 10,000 dumbells were considered at each nodal location in the model. The lid velocity is set to V = 1, the viscoelastic Weissenberg number We = 1 and Reynolds number of Re = 0.1. The simulation time of the movement is T = 1 in time increments of t = 0.0067 (N T = 150 snapshots).

The database consisted of the state vector (Eq. ( 23)) of the 100 nodes trajectories (excluding the node h = H , for which a no-slip condition v = 0 has been imposed). This database is split in 80 train trajectories and 20 test trajectories.

The net input and output size is N in = 5 and N out = 2N 2 in = 50. The number of hidden layers is N h = 5 with ReLU activation functions and linear in the last layer. It is initialized according to the Kaiming method [START_REF] He | Delving deep into rectifiers: surpassing human-level performance on imagenet classification[END_REF], with normal distribution and the optimizer used is Adam [START_REF] Diederik | Adam: a method for stochastic optimization[END_REF], with a weight decay of λ r = 10 The time evolution of the data L data and degeneracy L degen loss terms for each training epoch are shown in Fig. 9.

Results

Fig. 10 shows the time evolution of the state variables (position q, velocity v, internal energy e and conformation tensor shear component τ ) given by the solver and the neural net. There is a good agreement between both plots. Moreover, the proposed scheme is able to predict the time evolution of the flow for several snapshots beyond the training simulation time T = 1, as shown in the same figure.

Table 3 show the mean squared error of the data and degeneration loss terms for all the state variables of the Couette flow of an Oldroyd-B fluid. The results are computed separately as the mean over all the train and test trajectories using Eq. [START_REF] Grmela | Dynamics and thermodynamics of complex fluids. I. Development of a general formalism[END_REF]. Fig. 11 shows a box plot of the data error (MSE data ) for the train and test sets in the four studied architectures. The results of the structure-preserving neural network outperform the other two approaches even with noisy training data. The error of the unconstrained neural network is greater than one order of magnitude than our approach, proving the importance of the degeneracy conditions in the GENERIC formulation. Last, the naive black-box approach shows the worst performance of the four networks, as no physical restriction is considered. With respect to our previous work [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF], that employed a piece-wise linear regression approach, these examples show similar levels of accuracy, but a much greater level of robustness. For instance, this same example was included in the mentioned reference. However, in that case, the problem had to be solved with the help of a reduced order model with only six degrees of freedom, due to the computational burden of the approach. In our former approach, the GENERIC structure was identified by piece-wise linear regression for each of the few global modes of the approximation. So to speak, in that case, we learnt the characteristics of the flow. Here, on the contrary, the net is able to find an approximation for any velocity value at the 101 nodes of the mesh-say, fluid particles-without any difficulty. In this case, we are learning the

Conclusions

In this work we have presented a new methodology to ensure thermodynamic consistency in the deep learning of physical phenomena. In contrast to existing methods, this methodology does not need to know in advance any information related to balance equations or the precise form of the PDE governing the phenomena at hand. The method is constructed on top of the right thermodynamic principles that ensure the fulfillment of the energy dissipation and entropy production. It is valid, therefore, for conservative as well as dissipative systems, thus overcoming previous approaches in the field.

When compared with our previous works in the field (see Gonzalez et al. [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF]), the present methodology showed to be more robust, allowing us to find approximations for systems with orders of magnitude more degrees of freedom. This new approach is also less computationally demanding. For the double pendulum case, the snapshot optimization of the GENERIC matrices proposed in [START_REF] González | Thermodynamically consistent data-driven computational mechanics[END_REF] has a measured performance of 10 min per trajectory, which add up to 400 minutes considering the 40 studied trajectories, whereas our new neural-network approach trains in only 73.18 minutes. The computational time of the other examples is shown in Table 4.

The reported results show good agreement between the network output and the synthetic ground truth solution, even with moderate noisy data. We have also shown the importance of including the degeneracy conditions of the GENERIC formulation to the neural network constraints, as it ensures the thermodynamical consistency of the integrator. The structurepreserving neural network outperforms other naive black-box approaches, since the physical constraints act as an inductive bias, facilitating the learning process. However, the error can be reduced using several techniques:

• Database: As a general method of increasing the precision of an Euler integration scheme, the time step t can be decreased so the total number of snapshots is increased. On the contrary, the database will be larger, slowing the training process. The same way, the database can be enriched with a wider variety of cases, improving the net predictive capabilities.

• Integration Scheme: A higher order Runge-Kutta integration scheme could be introduced in Eq. (4) in order to get higher solution accuracy [START_REF] Wang | Runge-Kutta neural network for identification of dynamical systems in high accuracy[END_REF]. However, it requires several forward passes through the neural net for each time step, incrementing the complexity of the integration scheme and the training process. Additionally, GENERIC-based integration schemes have showed very good performance even for first-order approaches. [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF] • Net Architecture: To increase the computational power of the net, more and larger hidden layers N h can be added.

However, this could lead to a more over-fitted solution which limit the prediction power and versatility of the net. It also increases the computational cost of both the training process and the testing of the net.

• Training Hyperparameters: The neural networks trained in this work could be optimized using several hyperparameter tuning methods such as random search, Bayesian optimization or gradient-based optimization to get a more efficient solution.

Several open questions remain as a future work. A more exhaustive analysis can be performed to evaluate the influence of noisy data to the integrator evolution, in order to add robustness to the method and even predict wider simulation times using incremental learning [START_REF] Cybenko | Autoencoder-based incremental class learning without retraining on old data[END_REF][START_REF] Li | Learning without forgetting[END_REF].

CRediT authorship contribution statement

• Quercus Hernandez: software, investigation, validation • Alberto Badias: software, investigation, validation • David González: data, investigation, validation • Francisco Chinesta: methodology, validation • Elias Cueto: conceptualization, methodology, investigation, writing

Fig. 1 .

 1 Fig. 1. Representation of a single neuron (left) as a part of a fully connected neural net (right).

Fig. 2 .

 2 Fig. 2. Sketch of a structure-preserving neural network training algorithm.

Algorithm 1

 1 Pseudocode for the train algorithm. Load train database: z GT (train partition), t, L, M; Define network architecture: N in , N out = 2N 2 in , N h , σ j ; Define hyperparameters: η, λ d , λ r ; Initialize w i, j , b j ; for epoch ← 1, n epoch do for train_case ← 1, N train do Initialize state vector: z net 0 ← z GT 0 ;

Fig. 4 .

 4 Fig. 4. Loss evolution of data and degeneracy constraints for each epoch of the structure-preserving neural network training process of the double pendulum example.

  -5 and data loss weight of λ d = 10 2 . A multistep learning rate scheduler is used, starting in η = 10 -3 and decaying by a factor of γ = 0.1 in epochs 600 and 1200. The training process ends when a fixed number of epochs n epoch = 1800 is reached. The time evolution of the data L data and degeneracy L degen loss terms for each training epoch are shown in Fig. 4.

Fig. 5

 5 Fig.5. Time evolution of the state variables in a test trajectory of a double thermo-elastic pendulum using a time-stepping solver (Ground Truth, GT) and the proposed GENERIC integration scheme (SPNN). Since every variable has a vectorial character, both components are depicted and labeled as X and Y , respectively.

Fig. 8 .

 8 Fig. 8. Couette flow in an Oldroyd-B fluid.

Fig. 9 .

 9 Fig. 9. Loss evolution of data and degeneracy constraints for each epoch of the neural network training process of the Couette flow example.

  -5 and data loss weight of λ d = 10 3 . A multistep learning rate scheduler is used, starting in η = 10 -3 and decaying by a factor of γ = 0.1 in epochs 500 and 1000. The training process ends when a fixed number of epochs n epoch = 1500 is reached. The same parameters are considered also for the noisy database network (ν = 1%) and the unconstrained network.The black-box network training parameters are analogous to the structure-preserving network, except for the output size N out = N in = 5. Several network architectures were tested, and the lowest error is achieved with N h = 5 hidden layers and 25 neurons each layer.

Fig. 10 .

 10 Fig. 10. Time evolution of the state variables in five test nodes of a Couette flow using a solver (Ground Truth, GT) and the proposed GENERIC integration scheme (Net). The dotted vertical line represents the simulation time T = 1 of the training dataset.

Fig. 11 .

 11 Fig. 11. Box plots for the data integration mean squared error (MSE data ) of the Couette flow in both train and test cases.

  Pseudocode for the test algorithm.

	)
	Backward propagation;
	Optimizer step;
	end for
	Learning rate scheduler;
	end for
	Algorithm 2 Load test database: z GT (test partition), t, L, M;
	Load network parameters;
	for test_case ← 1, N test do
	Initialize state vector:

Table 1

 1 Mean squared error of the data loss (MSE

	State Variables		MSE data		MSE degen
			Train	Test	Train	Test
	q 1 [m]	X Y	1.95 • 10 -2 2 .72 • 10 -2	3.87 • 10 -2 8.21 • 10 -2	3.56 • 10 -8 4.74 • 10 -8	4.43 • 10 -8 5.77 • 10 -8
	q 2 [m]	X Y	2.04 • 10 -2 2 .72 • 10 -2	3.65 • 10 -2 3.73 • 10 -2	9.28 • 10 -8 3.55 • 10 -8	8.55 • 10 -8 5.01 • 10 -8
	p 1 [kg•m/s]	X Y	6.43 • 10 -4 1 .06 • 10 -3	1.33 • 10 -4 4.06 • 10 -3	4.00 • 10 -8 1.21 • 10 -7	7.08 • 10 -8 1.40 • 10 -7
	p 2 [kg•m/s]	X Y	4.88 • 10 -4 8 .47 • 10 -4	9.84 • 10 -4 1.79 • 10 -4	6.00 • 10 -8 9.76 • 10 -8	4.58 • 10 -8 1.20 • 10 -7
	s 1 [J/K]		1.21 • 10 -5	3.51 • 10 -5	1.31 • 10 -7	2.06 • 10 -7
	s 2 [J/K]		1.22 • 10 -5	3.18 • 10 -5	2.40 • 10 -7	2.95 • 10 -7

data ) and degeneracy loss (MSE degen ) for all the state variables of the double pendulum.

Table 2

 2 Mean squared error of the energy (MSE (E))

	and entropy (MSE (s)) of the double pendu-
	lum.
	Variable	Train	Test
	E [J]	7.99 • 10 -3	8.86 • 10 -3
	S [J/K]	6.52 • 10 -8	6.33 • 10 -8

Table 3

 3 Mean squared error of the data loss (MSE data ) and degeneracy loss (MSE degen ) for all the state variables of the Couette flow.

	State Variables	MSE data		MSE degen
			Train	Test	Train	Test
	q [-]	X	5.40 • 10 -6	6.29 • 10 -6	1.72 • 10 -7	1.96 • 10 -7
		Y	0 .00	0.00	0.00	0.00
	v [-]		3.23 • 10 -5	4.75 • 10 -5	1.19 • 10 -6	1.48 • 10 -6
	e [-]		7.85 • 10 -6	6.60 • 10 -6	7.06 • 10 -7	9.11 • 10 -7
	τ [-]		2.36 • 10 -5	1.26 • 10 -5	1.07 • 10 -6	1.31 • 10 -6

Table 4

 4 Computation training time of the proposed algorithm for the two reported examples in the noise free networks. It will be interesting, however, to study to what extent the employ of variational autoencoders, as in Bertalan et al.[START_REF] Bertalan | On learning Hamiltonian systems from data[END_REF], could help in solving more intricate models. Autoencoders help in determining the actual number of degrees of freedom needed to represent a given physical phenomenon.

	Example	Epoch Time	Total Time
	Double Pendulum	2.44 s/epoch	73.18 min
	Couette Flow	1.22 s/epoch	30.53 min
	behavior of fluid particles.

Acknowledgements

This project has been partially funded by the ESI Group through the ESI Chair at ENSAM Arts et Metiers Institute of Technology, and through the project 2019-0060 "Simulated Reality" at the University of Zaragoza. The support of the Spanish Ministry of Economy and Competitiveness through grant number CICYT-DPI2017-85139-C2-1-R and by the Regional Government of Aragon, grant T24_20R, and the European Social Fund, are also gratefully acknowledged.

✩ This project has been partially funded by the ESI Group through the ESI Chair at ENSAM ParisTech and through the project "Simulated Reality: an intelligence augmentation system based on hybrid twins and augmented reality" at the University of Zaragoza. The support of the Spanish Ministry of Economy and Competitiveness through grant number CICYT-DPI2017-85139-C2-1-R and by the Regional Government of Aragon, through the project T24_20R, and the European Social Fund, are also gratefully acknowledged.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.