Frédéric Besson
email: frederic.besson@inria.fr

Alexandre Dang
email: alexandre.dang@inria.fr

Thomas Jensen
email: thomas.jensen@irisa.fr

Securing Compilation Against Memory Probing

Keywords: Security and privacy → Formal security models, Software security engineering, Logic and verification, Secure Compilation, Side-channels, Dead Store Elimination, Register Allocation, Information-flow Preservation

A common security recommendation is to reduce the in-memory lifetime of secret values, in order to reduce the risk that an attacker can obtain secret data by probing memory. To mitigate this risk, secret values can be overwritten, at source level, after their last use. The problem we address here is how to ensure that a compiler preserve these mitigation efforts and thus that secret values are not easier to obtain at assembly level than at source level. We propose a formal definition of Information Flow Preserving program Transformations in which we model the information leak of a program using the notion of Attacker Knowledge. Program transformations are validated by relating the knowledge of the attacker before and after the transformation. We consider two classic compiler passes (Dead Store Elimination and Register Allocation) and show how to validate and, if needed, modify these transformations in order to be information flow preserving.

INTRODUCTION

Writing secure code is notoriously hard. It gets almost impossible if aggressive compiler optimisations may turn a secure source program into an insecure binary code. This can happen because of compiler bugs [START_REF] Yang | Finding and Understanding Bugs in C Compilers[END_REF] but also due to optimisations which exploit undefined behaviours of languages like C to generate code that does not match the expectation of the programmer [START_REF] Wang | Undefined Behavior: What Happened to My Code?[END_REF].

Using formal methods (e.g., abstract interpretation [START_REF] Blanchet | A Static Analyzer for Large Safety-critical Software[END_REF] or deductive verification [START_REF] Jacobs | VeriFast: A Powerful, Sound, Predictable, Fast Verifier for C and Java[END_REF]), it is possible to prove the absence of undefined behaviours of the source code. To avoid compiler bugs, verified compilers, e.g CompCert [START_REF] Leroy | Formal Certification of a Compiler Back-end or: Programming a Compiler with a Proof Assistant[END_REF] or CakeML [START_REF] Kiam Tan | A New Verified Compiler Backend for CakeML[END_REF], ensure that the observational behaviour of the source code is preserved by the compilation process.

However, security is a non-functional property that is not captured by the standard semantics of programs and therefore may not be preserved even if the program has a defined semantics and the compiler is correct. Recently, D'Silva et al. put focus on the "Correctness-Security Gap in Compiler optimizations" [START_REF] Vijay | The Correctness-Security Gap in Compiler Optimization[END_REF]. Similarly, the CERT at CMU recommends to "beware of compiler optimizations" 1 and provides the code example in Figure 1 to illustrate the problem. Here, the programmer limits the lifetime of a sensitive variable by overwriting its value using a memset. The semantics of As a result, a compliant C compiler is allowed to remove the call to memset and break the security of the code. The advisory also lists so-called compliant code with the caveat that the programmer [should] inspect the generated assembly code in the optimized release build to ensure that memory is actually cleared [...]. This is not very satisfactory.

Recently, Yang et al. [START_REF] Yang | Dead Store Elimination (Still) Considered Harmful[END_REF] have empirically studied the effectiveness of counter-measures and propose a generic memset which adapts to the platform to provide the maximum security. They also propose a secured version of Dead Store Elimination, the optimisation responsible for removing the memset. Simon et al. [START_REF] Simon | What you get is what you C: Controlling side effects in mainstream C compilers[END_REF] take another stance at the problem and propose to insert a compiler pass which secures the code by zeroing stack frames on context switches. These works provide in-depth empirical studies but do not formalise the end-to-end security property they intend to preserve or enforce.

In this work, we propose a formal definition of preservation of information leakage which, if verified by program transformations, ensures that the target code is not less secure than the source code, with respect to a passive but strong attacker able to read an arbitrary 1 CERT MSC06-C amount of memory. We stress that our purpose is not to enforce a security property; neither at the source level, nor at the target level. What we seek to identify is a general property that provides a formal account of how optimisations (such as the elimination of dead stores as in Figure 1) increase the information leakage and therefore render the optimised program insecure.

Our attacker model (ability to read arbitrary memory) is quite strong and one may wonder whether compiler optimisations preserve the information leakage. In this paper, we review some optimisations and show that some of them need to be adapted. Our contributions can be summarised as:

• We propose a notion of preservation of information leaks and assess its relevance against a list of simple transformations. • We show how to strengthen our initial proposal to capture information leaks due to partial observations. • We present sufficient conditions for the absence of information leaks which provide convenient reasoning principles. • We review two optimisations, dead store elimination and register allocation; and show how they can be adapted to preserve information leakage.

The rest of the paper is organised as follows. In Section 2, we propose a simple definition of preservation based on the notion of Attacker Knowledge. In Section 3, we present a strengthened definition which tracks finer information flows due to partial observations and is our definitive proposal of an Information Flow Preserving (IFP) transformation. Section 4 presents sufficient conditions for an IFP transformation which are easier to reason about in practice. Section 5 shows how the previous reasoning principles apply to two standard program optimisations: dead store elimination and register allocation. We discuss some extension of our formal model in Section 6. Related work is presented in Section 7 and Section 8 concludes.

INFORMATION-FLOW PRESERVATION

In this section, we present our formal definition of preservation of information leakage. We define our attacker model and propose a first version of what an Information Flow Preserving Transformation is. We evaluate this proposal against a series of program transformations, explain how it validates (or invalidates) simple transformations.

2.A Rationale for our Attacker Model

Our attacker model formalises the reason why reducing the lifetime of sensitive data is a good security practice. Basically, the risk to mitigate is that an attacker might be able to obtain sensitive data by reading the memory. For unsafe languages, e.g. C, any attacker code running in the same memory space can access arbitrary location and obtain sensitive data. Even if only trusted code is running, sensitive data can still be leaked through security vulnerabilities such as buffer overflows. In any case, to reduce the scope of these attacks, it makes sense to reduce the lifetime of sensitive data. Our attacker model is also relevant for type-safe languages, e.g. Java, where the memory is managed and reclaimed by a garbage collector. In this case, it is not possible for a source code program to read chunks of memory that is not accessible. Hence, if we consider an attacker being another type-safe program linked against the sensitive program, zeroing sensitive data would not be a relevant counter-measure. However, it still makes sense to protect against an attacker with physical access to the memory and therefore able to scan the memory content. In this context, it is worth overwriting sensitive data long before the occurrence of memory collection. Our attacker model covers both cases.

We consider an attacker which has total knowledge of the source and target code and can read arbitrary amount of memory. The attacker has a perfect knowledge of the memory layout and can target precisely the location of secret information. We will define what it means for a program transformation to be secure with respect to this attacker. Intuitively, we intend to ensure that an attacker does not get an advantage after the program transformation. In other words, the memory of the target program does not leak more information than the source code. In the following part, we formalise a simple but representative setting to illustrate these ideas.

2.B Formal setting

For simplicity, we consider deterministic, terminating executions and make the assumption that the attacker has only access to the final memory. As we explain in Section 6, these restrictions can be relaxed. In particular, it seems a reasonable trade-off for attackers to have memory access at the end of each function call. In order to model the information leak of an execution, we adapt the standard information-flow notion of Attacker Knowledge [START_REF] Askarov | Gradual Release: Unifying Declassification, Encryption and Key Release Policies (SP '07)[END_REF] to our simple setting. Definition 2.1 (Attacker Knowledge). Given a program p ∈ Program and a memory m ∈ Memory, the knowledge of the attacker (K (p, m)) is the set of initial configurations leading to the observation of the memory m. Formally, we have

K (p, m) = {c ∈ Config | (p, c) ⇓ m}
The notion of Attacker Knowledge characterises the uncertainty of the attacker about the initial configurations of the program. For a program p and a memory m, consider the following extreme cases. If K (p, m) = Config, the attacker has no information about the initial configuration because to her all of them compute the same memory m. On the contrary, if K (p, m) = {c} for some c, the attacker has perfect knowledge: she knows for sure that the program p is run from the configuration c. Hence, from a security standpoint, a bigger set K (p, m) is better. Definition 2.3 ensures that after a transformation the knowledge set of the attacker is getting bigger and therefore that she has less knowledge about the initial configurations leading to her observation. To simplify the definitions, we only consider source-to-source transformations. To capture transformation between source and target language, the previous definitions can be generalised in a straightforward way to distinguish between the source and target semantics.

2.B.2 Information Flow Preserving

(p, c) ⇓ m ∧ (p ′ , c) ⇓ m ′        ⇒ K (p, m) ⊆ K (p ′ , m ′ ) Based on

2.C Reality check

Definition 2.3 is elegant but it should also faithfully capture an intuitive notion of security i.e., reject program transformations deemed insecure and accept program transformations deemed secure. In the following, we consider the simple Program f of Figure 2 and several transformed versions of it. The Program f takes as input configuration the values of the variables x and y; zeroes the variable x and returns 0. As all the transformed programs compute the same value, i.e., 0, these are perfectly fine transformations from a compiler correctness point of view. In Figure 2, we mark with a ✓ the programs verifying Definition 2.2 and with a ✗ those which do not. The results are explained in details in Example 2.4.

Example 2.4. Consider Program g obtained by removing the assignment to the variable x from Program f. From a compiler point of view, this is a correct transformation because the value of x does not contribute to the result. As our motivating example shows (see Figure 1), removing dead store leaks more information to the attacker and should therefore be ruled out. Therefore we challenge our Definition 2.3 with the transformation from f to g which is insecure. Figure 3 Consider the initial configuration c = (x = 1, y = 0). The attacker observes the memory m = (x = 0, y = 0) for Program f and the memory m ′ = (x = 1, y = 0) for Program g. For Program f, because the value of the variable x is overwritten, the attacker has no information about its initial value. Therefore, the Attacker Knowledge is

K ( f , m) = {(x = 0, y = 0); (x = 1, y = 0)}
For Program g, the value of x is not overwritten and therefore the attacker has perfect knowledge of the initial configuration.

K (д, m ′ ) = {(x = 1, y = 0)}
As a result, for configuration c, we have

K ( f , m) ⊈ K (д, m ′ )
As a result, Definition 2.3 does not hold, therefore, as expected, the transformation from f to g is not information-flow preserving.

For the other examples of Figure 2, we give less formal explanation of whether they are information-flow preserving according to Definition 2.3. Consider Program h and Program i which both overwrite the value of variable x. Instead of the value 0, h uses the constant 1 and i uses the value of variable y. In terms of Attacker Knowledge, this difference vanishes as both assignments have the effect of erasing any information about the initial value of variable x. Furthermore, in both programs, the value of y is leaked, therefore, leaking it twice in Program i does not change the Attacker Knowledge. For any observed memory m, m ′ or m ′′ , we have

K ( f , m) = K (h, m ′ ) = K (i, m ′′ )
and the property holds for the Program h and Program i.

Program j introduces a direct information-flow from the initial value of variable x to variable a. It overwrites x but not a. As the variable a contains the exact same information as the initial variable of x, Program j is rejected. Program k also introduces the same direct information-flow but overwrites both copies i.e. a and x. In this case, the property holds. Program l has the same semantics as Program j but through an indirect information-flow instead of a direct information-flow. As the Attacker Knowledge is a semantics definition, it is immune to syntactic changes and therefore Program l is rejected for the exact same reason as Program j.

INFORMATION FLOW PRESERVATION

Definition 2.3 has the advantage of simplicity and matches our expectation for the programs of Figure 2. Nonetheless, as shown by the counter-example discussed in Section 3.A, it is still too permissive. In this section, we present (Definition 3.4), a refinement of Definition 2.3, which considers a hierarchy of attackers performing partial observations.

3.A Partial Observation and Information Leaks

Consider the programs f and f ' of Figure 4. Formally, the transformation verifies Definition 2.3. Yet, the transformed program f ' is intuitively easier to attack. In the following, we explain why an attacker with full observational power does not capture partial information leaks and show how to strengthen Definition 2. The intent of program f is best explained using a cryptographic analogy. The program f encrypts the message y using the one-time pad key x by performing an exclusive or with the message. Then, to make it hard for an attacker to recover the key, the program f overwrites the key x. The program f ' performs the same computation but does not overwrite the key. Intuitively, the transformation reduces security as the key x is leaked.

However, the transformation verifies Definition 2.3. Consider the scenario of Figure 4 where both programs run from the initial configuration c = (x = 1, y = 0). For Program f ' , the program computes the memory m ′ = (x = 1, y = 0, a = 1). Because the value of x is not overwritten, the attacker gains perfect knowledge of the initial configuration c: K ( f ′ , m ′ ) = {c}. For Program f, the value of x is erased and an attacker cannot directly observe its value. However, an attacker can still deduce the value of x from the values of a and y by solving the equation 1 = x ∧ 0. As the only solution is such that x = 1, we have: K ( f , m) = {c}. A similar reasoning can be made with any initial configuration and therefore the transformation satisfies Definition 2.3. Yet, this goes against our intuition that overwriting variable x makes the program more secure and should therefore be preserved by program transformations. In the following, we propose a strengthened security property which rules out the previous transformation.

3.B Discussion

In order to refine Definition 2.3, we first discuss why zeroing the key does not result in a less secure program. A first argument is that, for Program f ' , the key is present in-the-clear in the final memory whereas extracting it from Program f requires some computation. Following this path would require to equip the attacker with computational resources and relate the effort required to extract some piece of information. Another argument is that, for Program f, deducing 1 bit of key requires 2 bits of observation (both the plain-text and the cipher-text) whereas, for Program f ' , this can be done by observing directly 1 bit of the key. As a result, after the transformation, an attacker with a weaker observational power can reconstruct the key. In the following, we formalise this second approach which has the advantage of being purely logical and therefore model attackers without any computational limit.

3.C Partial Attacker Observation

We will represent attacker observations as the set of partial functions Addr ⇀ B ordered by the usual pointwise ordering:

o ⊑ o ′ iff ∀a ∈ dom(o), a ∈ dom(o ′ ) ∧ o(a) = o ′ (a)
We generalise this ordering to total functions. Hence, for an observation o and a memory m, o ⊑ m reads o is a partial observation of m. Definition 3.1 (Attacker Observation). Given a memory m and a number n of bits, we write obs (m, n) for the set of partial observations of the memory m limited to n bits. Formally, we have:

obs (m, n) = {o ∈ Addr ⇀ B | card (dom(o)) = n ∧ o ⊑ m}
K m n (p, o) = o ⊑m ′ K (p, m ′ )
If the observation is total (n = card (Addr)), we get the obs (m, n) = {m} and the only memory such that m ⊑ m ′ is m itself. Under those conditions, Definition 3.2 coincides with Definition 2.1.

K m n (p, o) = K (p, m)
In the general case, as the observation o is partial, all the memories m such that o ⊑ m are indistinguishable to the attacker and are potential results of the program. As the result, the Attacker Knowledge obtained from the observation o is the union of the knowledge obtained for each potential memory m.

3.D Hierarchy of Attackers

Using Definition 3.1, we define a stronger notion of informationflow preservation.

Definition 3.3. Given an attacker observing n bits of memory, a program p is less secure than a program p ′ (written p ≼ n p ′ ) if their respective Attacker Knowledge are ordered such that ∀c, m, m ′ we have: Using the previous ordering of programs, we are ready to define our notion of information-flow preservation. As before, consider that both programs are run from the initial configuration c = (x = 1, y = 1). In f ' , the attacker performs a 1-bit observation and picks the value of x i.e., the key that is not erased. From this observation, she gets perfect knowledge of the value of x but obviously no information about y.

(p, c) ⇓ m ∧ (p ′ , c) ⇓ m ′        ⇒ ∀o ′ .∃o.K m n (p, o) ⊆ K m ′ n (p ′ , o ′ ) Definition 3.
K m ′ 1 ( f ′ , (x = 1)) = {(x = 1, y = 0); (x = 1, y = 1)}.
To prove that the transformation is secure, Definition 3.4 mandates that there should exist a stronger 1-bit observation for program f. As shown by Figure 5, this is not possible. We examine in turn each of the possible 1-bit observations.

• As the variable x is overwritten, no knowledge at all is derived from its observation

K m 1 ( f , (x = 0)) = Config. • Observing y gives no knowledge over x K m 1 ( f , (y = 0)) = {(x = 0, y = 0), (x = 1, y = 0)}.
• Observing a gives some information about the possible initial configurations. Yet, the knowledge set is not included in

K m ′ 1 ( f ′ , (x = 1)). K m 1 ( f , (a = 1)) = {(x = 1, y = 0), (x = 0, y = 1)}.
As a result, the transformation is not Information Flow Preserving.

SUFFICIENT CONDITIONS FOR IFP

The next step is to identify a set of sufficient conditions to prove that a program transformation is IFP (Section 4.A and Section 4.B) and use them to study two compiler optimisations: Dead Store Elimination (Section 5.A) and Register Allocation (Section 5.B).

4.A Equivalent Programs

As our IFP property is an end-to-end property, it is independent of syntactic differences between programs. As a result, if two programs always compute the same memories, the transformation is IFP. Definition 4.1 (Equivalent programs). Two programs p and p ′ are equivalent (written p ≡ p ′ ) if for every configuration c ∈ Config they compute the same memory. Formally, we have:

p ≡ p ′ ⇔ ∀c, m, m ′ .((p, c) ⇓ m) ∧ ((p ′ , c) ⇓ m ′ ) ⇒ m = m ′
(p, c) ⇓ m ∧ (T (p), c) ⇓ m ′        ⇒ ∀o ′ .∃o.K m n (p, o) ⊆ K m ′ n (T (p), o ′ )
First, notice that for two equivalent programs, the Knowledge of the attacker is the same, i.e., we have

K m n (p, o ′ ) = K m n (T (p), o ′ ). (1)
For every observation o ′ , we take o to be o ′ . It then remains to prove that

K m n (p, o ′ ) ⊆ K m ′ n (T (p), o ′ ).
As p ≡ T (p), we have m = m ′ and by Equation 1 the property holds. □

4.B Constant and Matching Addresses

Definition 3.3 relates program executions using the general notion of Attacker Knowledge. In practice, it is often more convenient to characterise explicitly the relation between the produced memories before and after an IFP transformation.

Note that a memory address whose value is always constant contains no knowledge about the input configuration. This remark is useful to restrict the relation between memories to only relevant (i.e., non constant) addresses. Definition 4.3 (Constant address). For a transformed program p ′ , the set of constant addresses is the set of addresses which return the same value for every execution of the program. Formally, we have:

cst (p ′ ) = a ∃v.∀c.(p ′ , c) ⇓ m ⇒ m(a) = v
A second remark is that if two addresses a and a ′ always contain the same value, their knowledge is identical and therefore an observation of a ′ can be matched with an observation of a. Definition 4.4 (Matching Address). Given two programs p and p ′ , the addresses a and a ′ are in matching correspondence (written p a ≡ a ′p ′ ) if for any execution, they contain the same value. Formally, we have:

p a ≡ a ′p ′ ⇔ ∀c.(p, c) ⇓ m ∧ (p ′ , c) ⇓ m ′ ⇒ m ′ (a ′ ) = m(a)
By combining Definition 4.3 and Definition 4.4, we obtain a sufficient condition for an IFP transformation. Intuitively, a transformation is IFP if for every address a ′ that is not constant for p ′ there exists a matching address a in p. This is stated formally in Theorem 4.5.

Theorem 4.5. The transformation from program p to program p ′ is IFP if

∀a ′ .a ′ cst (p ′ ) ⇒ ∃a.p a ≡ a ′p ′ ⇒ ∀n.p ≼ n p ′ .
Proof. A comprehensive proof can be found in Appendix A. The following proof sketch presents the main arguments. By expressing our theorem in terms of attacker knowledge our goal is to prove the following:

∀n.∀o ′ .∃o.K m n (p, o) ⊆ K m ′ n (p ′ , o ′ )
We proceed by induction on n:

• n = 0. This part is straightforward since the attacker can not observe any bits of the returned memory. The only observation that can be made by the attacker is ⊥, the observation with an empty domain. Hence the attacker knowledge we obtain is the set of all the possible inputs Config.

∀ p m, K m 0 (p, ⊥) = Config With this relation we can easily prove our goal for n = 0.

• n → n + 1. We slightly simplify the notation

o ⊑ m K (p, m) → o K (p, m)
for readability. The relation o ⊑ m is implicit with the memory used as the second argument of attacker knowledge.

A relation between attacker knowledge for n and n + 1 exists and is the following (more details in Appendix A): Starting from our goal proof we develop the following proposition with the help of our induction hypothesis:

K m[a→v] n+1 (p, o[a → v]) = K m[a→v] n (p, o) \ o K (p, m[a ↛ v]) We define m[a → v]
K m[a→v] n+1 (p, o[a → v]) ⊆ K m ′ [a ′ →v ′ ] n+1 (p ′ , o ′ [a ′ → v ′ ]) ⇔ K m[a→v] n (p, o) \ o K (p, m[a ↛ v]) ⊆ K m ′ [a ′ →v ′ ] n (p ′ , o ′ ) \ o ′ K (p ′ , m′ [a ′ ↛ v ′ ]) ⇔ o K (p, m[a ↛ v]) ⊇ o ′ K (p ′ , m′ [a ′ ↛ v ′ ])
Recall that in our theorem we have the hypothesis that a ′ is either a constant or a matching address. We prove the proposition in both cases:

• a ′ ∈ cst(p ′ ). In this case we know that v ′ is a constant such that for any input the returned value at address a ′ is always v ′ . This implies that:

o ′ K (p ′ , m′ [a ′ ↛ v ′ ]) = ∅
Indeed, since for all executions the value at a ′ is always v ′ , the memories referred by m′ [a ′ ↛ v ′ ] do not exist. With this equality we can easily prove our proposition.

• p a ≡ a ′ p ′ . We know that a ′ is a matching address to a, therefore for every configuration of p and p ′ leading to m and m ′ we have m(a ′ ) = m(a) = v. We use the definition of attacker knowledge on our proposition which gives:

{c | (p, c) ⇓ m[a ↛ v]} ⊇ {c | (p ′ , c) ⇓ m′ [a ′ ↛ v]} (2)
We need to prove that for every configuration c where (p ′ , c) ⇓ m′ with m′ (a ′ ) v then we also have (p, c) ⇓ m with m(a) v. This is straightforward since we know that m(a) = m′ (a ′ ) then having m′ (a ′ ) v necessarily implies m(a) v. We have proved our proposition and Theorem 4.5. □ Theorem 4.5 provide a sound characterisation of an IFP transformation. However, it is not complete. This is illustrated by the following IFP transformation x = y → x =∼ y which does not meet the pre-conditions of Theorem 4.5. Yet, the pre-conditions are sufficiently general to cover sophisticated program transformations such as Register Allocation (see Section 5.B).

IFP TRANSFORMATIONS

In this section we focus on two compilation steps which are not a priori IFP: Dead Store Elimination (DSE) and Register Allocation (RA). We show how to make these two transformations IFP and we prove them sound using two different proof techniques. First we present a modified DSE and prove directly that its algorithm is sound. Second, we show that RA can be made IFP by resorting to an a posteriori approach using a certified validator which checks the IFP condition for a particular program and its RA transformation.

5.A Dead Store Elimination

DSE is the archetypical program transformation that is not IFP. Informally, a dead store is a memory write which is provably unnecessary to compute the program result. Hence, from an optimisation point of view, it is a perfectly legal transformation to remove a dead store instruction. Security-wise, as shown by our motivating example (see Figure 1) and Program g of Figure 2, the transformation does not preserve information-flows. Indeed, an attacker may gain more knowledge by observing the value that is not overwritten due to the removal of a dead store.

A drastic solution would be to disable this optimisation. We propose a modified Dead Store Elimination optimisation based on a revised and strengthened notion of dead store.

5.A.1 Liveness based DSE.

Dead stores are typically identified using a liveness analysis [15, 2.1.4]. A liveness analysis is a classic backward program analysis. For each program point, it computes an over-approximation of the live variables i.e. variables that are necessary to compute the program result. Dually, dead variables are those variables that are not live and a dead store is a memory write (x = e) where the variable x is dead. Therefore, a classic DSE performs a liveness analysis and removes dead stores.

Our abstract model of programs does not feature program points. Instead we use the end-to-end property of a liveness analysis: Definition 5.1 (Sound Liveness). Given a program p, an end-to-end liveness analysis is a pair of sets of addresses (V ,W ) ⊆ Addr ×Addr such that

∀c, c ′ , m, m ′ . (p, c) ⇓ m ∧ (p, c ′ ) ⇓ m ′        ⇒ c ≃ V c ′ ⇒ m ≃ W m ′ where e ≃ V e ′ ⇔ ∀x ∈ V .e (x ) = e ′ (x )
for (e, e ′ ) ∈ Config × Config or (e, e ′ ) ∈ Memory × Memory.

Definition 5.1 says that if two input configurations c and c ′ agree on the values of the input live addresses V then the output memories agree on the output live addresses W . For the purpose of optimisation, the set W contains a minimum set of addresses i.e., only those needed by the return statement of the program; the set V being computed by a backward fixpoint iteration (see [15, 2.1.4] for details).

As DSE removes the dead statements, it keeps unchanged the values of live addresses. As a result, in our formal model, a DSE transformation can be characterised by Definition 5.2. Definition 5.2 (Dead Store Elimination). Given a program p and a liveness analysis (V ,W ) of p. The DSE transformation from p to p ′ is correct if p ′ is indistinguishable from p for all the live variables W . Formally, we have:

∀c.(p, c) ⇓ m ∧ (p ′ , c) ⇓ m ′ ⇒ m ≃ W m ′
This specification is partial. In particular, it says nothing about the dead variables but this is enough to conclude that these dead variables may leak information and violate our IFP property.

5.A.2 Shadow Store Elimination.

In order to get an IFP DSE, we propose to exploit Theorem 4.2 and therefore have a dead store elimination producing an equivalent program according to Definition 4.1. Interestingly, this result can be obtained in a non-intrusive way by only slightly modifying the initial condition of the liveness analysis. Indeed, it is sufficient to impose that, at the end of the program, every address is live. As liveness analysis is backward, this can always be done. The effect of this modification is that only dead stores that are shadowed by a following store at the same address can be safely removed. Said otherwise, for a given address, the last store needs to be retained (whether it is dead or not). However, the penultimate dead store can be safely removed as its value is overwritten later on.

5.B Register allocation

Register Allocation (RA) is an essential compiler pass which makes sure that the target program abides to the constraints of the underlying architecture. In particular, RA compiles source programs using an unlimited number of variables into target programs using a limited number of hardware registers. RA is based on sophisticated algorithms and heuristics e.g., Iterated Register Coalescing [START_REF] George | Iterated Register Coalescing[END_REF], which aim at optimising the usage of registers. As a last resort, when no more registers are available, the RA algorithm performs a socalled register spill i.e., the value of the register is temporarily saved in memory so that the register can be reused for other purpose. The symmetric operation is register reloading consisting in restoring the saved value.

Example 5.4. Figure 6 illustrates the working of a RA algorithm for an architecture with only two registers. The source code on the left computes the bitwise xor of text with some salt and some key. Note that the source code is implicitly erasing the key by reusing the variable c for computing the program result. The target code on the right is the result of a typical RA algorithm. The architecture has only two registers i.e., r1 and r2. Moreover, the architecture imposes that a xor operation can only be performed with register operands.

Initially, the allocator maps the variable c to the register r2 and the variable salt to the register r1. However, to perform a xor with text , the allocator needs to reuse a register. To do that, the register r2 is spilled in memory in the variable c so that r2 can be reused for storing the program argument text . For the second xor, the variable c is reloaded in register r1.

From a security point of view, the register allocation of Figure 6 is not an IFP transformation. The culprit is the spilling of register r2. Whereas the value of r2 is eventually overwritten, its copy remains in memory in variable c. Compared to DSE, the introduction of the information flow leak is not due to an aggressive optimisation but is really needed for a successful RA.

To ensure and enforce that RA is an IFP transformation we propose a variation of the translation validation approach [START_REF] Pnueli | Translation Validation[END_REF].

5.B.1 Translation

Validation of Register Allocation. The design of our validator is guided by Theorem 4.5 which states sufficient conditions for an IFP transformation. Definition 5.5 states the soundness conditions for such a validator. Definition 5.5 (Sound Validator). Let p be a program and p ′ be a transformed program. A sound validator is such that IFPok(p, p ′ ) = true entails that

∀a ′ .a ′ cst (p ′ ) ⇒ ∃a.p a ≡ a ′p ′
In other words if a ′ is not a constant address of p ′ , then there must exist a matching address a in p. As the conditions quantify over all addresses, a naive approach i.e., enumerating over every address, is not feasible. To get an effective solution, notice that RA keeps the memory mostly unchanged with the exception of a limited number of locations which are program variables, registers and spilled locations. For our formal model, consider that the set of addresses is partitioned into a set of named addresses NAddr, affected by register allocation, and unamed addresses UAddr untouched by register allocation. As a result, the validator needs to establish a precise mapping between named addresses but only needs to check that other addresses have the exact same value before and after the transformation.

The validator performs a constant analysis of the transformed program thus identifying addresses which do not leak any information. Constant analysis is a standard program analysis. For our context, Definition 5.6 states the soundness guarantee of such a constant analysis. Definition 5.6 (Soundness of Constant Analysis). For a program p, a sound constant analysis cst ♯ (p) ∈ NAddr ⇀ B is a partial function mapping named addresses to a constant values such that

dom(cst ♯ (p)) ⊆ cst (p)
A matching analysis is less standard but existing validators for RA [START_REF] Huang | Catching and Identifying Bugs in Register Allocation[END_REF][START_REF] Rideau | Validating Register Allocation and Spilling[END_REF] implement a similar technique. Definition 5.7 states the soundness requirements. Definition 5.7 (Sound Matching Analysis). For programs p and p ′ , a sound matching analysis

match ♯ ∈ Program × Program ⇀ (NAddr ⇀ NAddr)
is such that if match ♯ (p, p ′ ) is defined we have that • all the unnamed addresses have the same value in both p and p ′ (∀a ′ ∈ UAddr, p a ′ ≡ a ′p ′ ); • and the obtained mapping φ ∈ NAddr ⇀ NAddr relates matching addresses (∀a ′ ∈ NAddr, p φ (a ′ ) ≡ a ′p ′ ).

Using the previous analyses, we can construct a sound IFP validator for RA. Theorem 5.8 (Sound Validator). Let p be a program and p ′ obtained by an untrusted RA. Suppose a RA validator defined by

IFPok(p, p ′ ) = ∃φ.match ♯ (p, p ′ ) = φ ∧ ∀a ′ cst ♯ (p ′ ).∃a.φ(a ′ ) = a If IFPok(p, p ′ ) = true then the transformation is IFP.
Proof. By Theorem 4.5, it amounts to proving

∀a ′ .a ′ cst (p ′ ) ⇒ ∃a.p a ≡ a ′p ′
By hypothesis, the matching analysis succeeds such that match ♯ (p, p ′ ) = φ for some φ. The proof is case analysis over the address a ′ .

• If a ′ ∈ UAddr, by Definition 5.7 and because the validator succeeds, we have p a ′ ≡ a ′p ′ and therefore the property holds.

• If a ′ ∈ NAddr ∧ a ′ ∈ dom(cst ♯ (p ′ ))
. By Definition 5.6, we have a ′ ∈ cst (p ′ ) and, by contradiction, the property holds. Example 5.9. The matching analysis for the programs of Figure 6 is shown Figure 7. For clarity, variable names are prefixed by a S_ in the source code and by a T_ in the target code. A matching is attached to a pair of program points and we write φ i, j for a matching φ where i is a program point of the source program and j is a program point of the target program. As both programs execute, the matching is updated. In Figure 7, a dashed line links an update of the mapping with the instruction responsible for the update.

• If a ′ ∈ NAddr ∧ a ′ cst ♯ (p ′ ),
At the start of both programs the parameters are stored at the same address. As the result, we have the initial matching φ 1,1 = [T _text → S_text]. After the first two instructions, we have the matching φ 3,3 = [r 2 → S_C, r 1 → S_salt,T _text → S_text]. The next two instructions of the target code correspond to the spilling of register r 2 and loading of variable T _text and have no counterpart in the source code. Therefore, we get the matching φ 3,5 = [r 2 → S_text, r 1 → S_salt,T _c → S_c,T _text → S_text]. At the entry of the first xor, the operator and the operands do match and after the instruction the matching is updated and we get:

φ 4,6 = [r 2 → tmp, r 1 → S_salt,T _c → S_c,T _text → S_text]
At the instruction 8 of the transformed program, the register r 1 is mapped to S_C. However, as the source variable S_c is overwritten, the variable T _c has no more mapping with a source variable hence we remove it from φ. We get the final mapping

φ 5,9 = [r 2 → tmp, r 1 → S_c,T _text → S_text]
At the end of program, the variable T _c has no mapping. This is the witness of the information flow leak and the reason why the validator rejects the transformation.

5.B.2 From Translation

Validation to IFP RA. As register allocation is, in general, not IFP, there is very little chance that our validator would ever succeed. However, as illustrated by Example 5.9, the missing mappings of the validator indicate the information flow leaks. Closing them can be done by zeroing all those variables or registers at the end of the program. For the program of Figure 7, this amounts to adding the assignment T_c = 0. After the transformation, the matching analysis would be unchanged. However, the constant analysis would detect that the value of T _c → 0 i.e., leaks no information, and therefore prove that the transformation is IFP.

This approach is non-intrusive and compatible with any RA algorithm. Patching the information leaks a posteriori may not be the most efficient approach and an IFP aware RA algorithm could make a wiser reuse of spilling locations. Yet, practical experiments are needed to assess the effect on the quality of the generated code. In the worst case, all spilled locations and registers would need to be overwritten. In practice, as calling conventions mandate calleesaved registers to be restored this worst-case is fortunately never reached.

DISCUSSION

For the sake of clarity, our formal model is simple. Yet, the result of this paper could be generalised to more complex contexts e.g., where programs executions are not necessarily terminating and the attacker is not limited to a single observation at the end of the program.

For Termination, Definition 3.3 takes as hypothesis that both programs terminate and is therefore termination sensitive. If a nonterminating program p is transformed into a terminating program p ′ , the property is vacuously true whereas the transformed program potentially leaks more information. To deal with this issue, it is sufficient to require that p ′ only diverges if p already does. As most compiler optimisations preserve termination, this is a mild requirement which always holds in practice.

For the Attacker model, our assessment is that an omniscient Attacker able to observe memory at any time would preclude too many optimisations. A more practical scenario would be to define pre-determined observation points, e.g., function exits, which can be matched between a program and its transformation. In this case, our IFP definition would need to be generalised. In particular, if the attacker possesses a memory, each observation would monotonically increase its knowledge.

Even after refinement, our attacker model is still quite strong and does not model how difficult it is to precisely probe the memory. For instance, it considers equally insecure a program leaking a key spread over the whole memory or the same key stored in adjacent memory locations. Similarly, for our attacker model, duplicating (secret) information does not increase the information flow leak whereas, intuitively, it makes practical attacks more feasible. An interesting extension would be to equip the attacker with probability distributions modelling how hard it is to successfully and precisely read a given memory location.

RELATED WORK

The problem of secure erasure of secrets has been studied from an information-flow perspective. Chong and Myers [START_REF] Chong | Language-Based Information Erasure[END_REF] propose semantics foundations for defining erasure policies. A main insight behind that work is that erasure can be seen as the dual of declassification [START_REF] Sabelfeld | Declassification: Dimensions and principles[END_REF]. Later on, Chong and Myers [START_REF] Chong | End-to-End Enforcement of Erasure and Declassification[END_REF], but also Hunt and Sands [START_REF] Hunt | Just Forget It: The Semantics and Enforcement of Information Erasure[END_REF], propose type-systems for verifying erasure policies of the source code. Askarov et al. enforce the erasure policies in the presence of so-called write-once locations which cannot be overwritten. The key insight is to store encrypted data in write-once locations and simulate erasure by the erasure of the cryptographic key. In our work, we do not consider erasure policies but ensure that program transformations do not introduce information leaks. Our model of attacker is currently too idealised to ensure that IFP transformations preserve the rich erasure policies defined by the aforementioned works. An attacker observing all the intermediate memory states would probably suffice but would preclude too many transformations.

Recently, D'Silva et al. [START_REF] Vijay | The Correctness-Security Gap in Compiler Optimization[END_REF] drew the attention to the correctnesssecurity gap, explicitly mentioning dead-store elimination as a problematic transformation. Deng and Namjoshi [START_REF] Deng | Securing a Compiler Transformation[END_REF] ensure that a modified DSE transformation preserves their information flow notion of leaky triple. This notion does not capture the amount of information leak. For instance, it ensures that a non-interferent program is transformed into another non-interferent program but allows a leaky program to be transformed into a more leaky program. Our IFP property, based on the notion of Attacker Knowledge, is stronger and rules out this possibility: it ensures that the transformed program is no more leaky. Yang et al. [START_REF] Yang | Dead Store Elimination (Still) Considered Harmful[END_REF] survey the workarounds used by developers to prevent compilers to remove security-sensitive dead-stores. They also propose a secure DSE implementation for LLVM which verifies our IFP property. Simon, Chisnall and Anderson [START_REF] Simon | What you get is what you C: Controlling side effects in mainstream C compilers[END_REF] investigate how compilers may break the security of cryptographic code. They propose compiler support for the secure erasure of secret which improves the security of the generated code.

Recently, Barthe, Grégoire and Laporte [START_REF] Barthe | Secure compilation of side-channel countermeasures: the case of cryptographic "constant-time[END_REF] have proposed a general framework to reason about the preservation of information leakage. They use an instrumented semantics which explicitly leaks information to the attacker and propose reasoning principles based on so-called 2-simulations. It is unclear whether our attacker model fits their framework because our attacker performs an arbitrary, non-deterministic, observation of the memory. Our reasoning principles are simpler because our necessary conditions for IFP (see Theorem 4.5) can be enforced using traditional simulation proof techniques for compiler correctness [START_REF] Leroy | Formal verification of a realistic compiler[END_REF].

CONCLUSION

Compiler optimisations introduce information leaks which make them unreliable for security crtitical code. In this work, we propose a formal definition of Information Flow Preserving (IFP) transformation using an attacker with the ability to observe memory. A transformation is IFP if an attacker does not gain any advantage by observing the final memory after running the transformed program. We have shown that classic optimisations such as Dead Store Elimination (DSE) and Register Allocation (RA) are in general not IFP transformation. To our knowledge, this is the first time that DSE is shown insecure based on an Information Flow argument. The good news is that, for DSE and RA, modifying the transformations to be IFP is feasible in a rather non-intrusive way.

In theory, optimisations that improve the instruction sequence but preserve memory writes (e.g. constant propagation or common expression elimination) are IFP transformations. As future work, we intend to investigate more optimisations and review their actual implementation. On a more practical side, we intend to evaluate the difficulty of modifying a realistic compiler to make it IFP but also the loss of efficiency this modification may incur. Another limitation to address is that we only consider whole-program optimisations. Further work should examine to what extent it is possible to obtain a more compositional model that can reason about individual functions and in which an attacker can observe memory in-between function calls.

A PROOF OF THEOREM 4.5

We imagine an execution of p and p ′ resulting in memories m and m ′ respectively. We keep for later use the hypothesis on p and p ′ : ∀a ′ .a ′ cst (p ′ ) ⇒ ∃a.p a ≡ a ′p ′ (H)

Our goal is to prove: We do an induction on n ∈ N.

• n = 0. This part is straightforward since the attacker can not observe any bits of the returned memory. The only observation that can be made by the attacker is ⊥, the observation with an empty domain. Hence the attacker knowledge we obtain is the set of all the possible inputs Config.

∀ p m, K m 0 (p, ⊥) = Config
With this relation we can easily prove our goal for n = 0.

• n → n + 1. The remaining part to prove for our induction is that if our theorem is true for n it is also true for n + 1. Before describing the hypothesis and goals we search for a relation between attacker knowledge observing n and n + 1 bits.

We use the notation o[a → v] for an observation of n + 1 bits composed of the observation of n bits o and the observation of one bit at address a which returns v.

We detail the definition of an attacker knowledge of n bits on the memory m[a → v] where m(a) = v. The attacker knowledge is built with the observation o where we assume a dom(o). We slightly simplify the notation from From (4) to [START_REF] Chong | End-to-End Enforcement of Erasure and Declassification[END_REF] we only apply the definition of attacker knowledge for an observation o[a → v] of n + 1 bits.

With this new relation we can prove that if the property is true for n then it is also true for n + 1. We develop the following proposition with the help of the induction hypothesis (H):

K m[a→v] n+1 (p, o[a → v]) ⊆ K m ′ [a ′ →v ′ ] n+1 (p ′ , o ′ [a ′ → v ′ ]) ⇔ K m[a→v] n (p, o) \ o K (p, m[a ↛ v]) ⊆ K m ′ [a ′ →v ′ ] n (p ′ , o ′ ) \ o ′ K (p ′ , m′ [a ′ ↛ v ′ ]) ⇔ o K (p, m[a ↛ v]) ⊇ o ′ K (p ′ , m′ [a ′ ↛ v ′ ])
We remind that in our theorem we have the hypothesis that a ′ is either a constant or a matching address. We prove the proposition in both cases:

• a ′ ∈ cst(p ′ ). In this case we know that v ′ is a constant such that for any input the returned value at address a ′ is always v ′ . Meaning that:

o ′ K (p ′ , m′ [a ′ ↛ v ′ ]) = ∅
Indeed since for all executions the value at a ′ is always v ′ , the memories referred by m′ [a ′ ↛ v ′ ] do not exist. With this equality we can easily prove our proposition.

• p a ≡ a ′ p ′ . We know that a ′ is a matching address to a therefore for every configuration of p and p ′ leading to m and m ′ we have m(a ′ ) = m(a) = v. We use the definition of attacker knowledge on our proposition which gives:

{c | (p, c) ⇓ m[a ↛ v]} ⊇ {c | (p ′ , c) ⇓ m′ [a ′ ↛ v]} (3)
We need to prove that for every configuration c where (p ′ , c) ⇓ m′ with m′ (a ′ ) v then we also have (p, c) ⇓ m with m(a) v. This is straightforward since we know that m(a) = m′ (a ′ ) then having m′ (a ′ ) v necessarily implies m(a) v. We have proved our goal and Theorem 4.5.

1 4 / 5 } 6 memsetFigure 1 :

 4561 Figure 1: Non-compliant erasure

2.B. 1

 1 Knowledge of the Attacker. We consider a finite, bit-addressable memory Memory = Addr → B where B = {0, 1} is the set of boolean values. A program p ∈ Program runs from an initial configuration c ∈ Config and produces a final memory m ∈ Memory. This is written (p, c) ⇓ m. For the sake of the presentation, we only consider terminating programs i.e., ∀c.∃m.(p, c) ⇓ m. We discuss in Section 6 how to lift this restriction.

  exhibits this example.

Figure 3 :

 3 Figure 3: Counter-example to Knowledge inclusion

  memory c = (x = 1, y = 0)

Figure 4 :

 4 Figure 4: Paradox: not erasing encryption keys preserves information leakage

Definition 3 .

 3 2 generalises the notion of Attacker Knowledge for an observation. Definition 3.2 (Attacker Knowledge from Observation). For a given number of bits n and a memory m, the attacker knowledge for a given observation o ∈ obs (n, m) is defined by

  3 requires an alternation of quantifiers to model the fact that the Attacker of program p ′ chooses the partial observation o ′ it intends to perform and challenges the Attacker of program p to perform a more powerful observation o with a smaller knowledge set than the observation o ′ . Note also that if the observations are total (n = card (Addr )), Definition 3.3 coincides with Definition 2.2 because we necessarily have o ′ = m ′ and o = m.

Definition 3 . 4 (Figure 5 :

 345 Figure 5: not IFP for an attacker with 1 bit of observation

Theorem 4 .

 4 2 states that a transformation mapping a program to an equivalent program is IFP. Theorem 4.2. Let T : Program → Program. We have (∀p, p ≡ T (p)) ⇒ IFP (T ) Proof. By Definition 3.3 and Definition 3.4, we have to prove that:

  as a memory satisfying m(a) = v. The observation o has domain of size n and a dom(o). o[a → v] is the composition of o and the function a → v hence its domain has cardinal n + 1. Therefore attacker knowledge of n + 1 on m[a → v] with observation o[a → v] is equal to the attacker knowledge of n with observation o minus the union of all the possible configurations leading to a potential memory m satisfying both o ⊑ m and m(a) v.

Theorem 5 . 3 .

 53 Given a program p and a liveness analysis (V ,W ) of p such that W = Addr, DSE produces a program p ′ that is IFP. Proof. By Theorem 4.2, it suffices to prove p ≡ p ′ . Suppose that (p, c) ⇓ m and (p ′ , c) ⇓ m ′ . By Definition 5.2 of a DSE transformation, we have m ≃ Addr m ′ i.e., m = m ′ . As a result, Theorem 5.3 holds.□

e l o a d c i n r 1 8 r 1 = r 1 ^r 2 9 r e t u r n r 1 Figure 6 :

 81916 Figure 6: Original program (left) After register allocation (right)

Figure 7 :

 7 Figure 7: Matching analysis example

  ∀n.p≼ n p ′ ⇔ ∀n.∀o.∃o ′ . K m n (p, o) ⊆ K m ′ n (p ′ , o ′ )

  o ⊑ m K (p, m) to o K (p, m) for readability. The relation o ⊑ m is implicit with the memory used as the second argument of attacker knowledge.

  o[a → v]) (5)Equation (1) is simply the definition of Attacker Knowledge. m is the possible memories which can have o as a partial function: o ⊑ m.From (1) to (2) we split Attacker Knowledge into two. First the configurations whose execution results in memory m where the value at a is v ( m[a → v]). Second the memories where the value ata is different from v ( m[a ↛ v]).To go from (2) to (3) we pass the term o K (p, m[a ↛ v]) to the left side of the equality turning the operator ∪ into \. This is true because o K (p, m[a ↛ v]) and o K (p, m[a → v]) are disjoint sets. Indeed since our programs are deterministic the execution of a configuration c can not lead to both m[a → v] and m[a ↛ v]. From (3) to (4) we use the equality: o K (p, m[a → v]) = o K (p, m) Instead of considering memories m[a → v] that satisfy both o ⊑ m and m(a) = v, we only consider memories that satisfy o[a → v] ⊑ m. For reminder, o[a → v] is the composition of o and the function that associates a to the value v. Therefore the equality is true since memories satisfying o[a → v] ⊑ m also satisfy the two conditions on m[a → v].

  Transformation as Increase of Attacker Knowledge . Intuitively, a program transformation is secure against our attacker if her knowledge does not increase by examining the memory after the program transformation. Definition 2.2 formally defines a partial order on programs such that p ≼ p ′ is read as the program p ′ leaks no more information than For convenience, we say that p leaks more information than p ′ , keeping in mind that the relation is reflexive.

	Original program
	f : x=0; return 0;
	Dead store elimination
	g: skip ; return 0;	✗
	Overwrite with constant value
	h: x=1; return 0;	✓
	Overwrite with leaked value
	i : x=y; return 0;	✓
	Leaking through direct flow
	j : a=x; x=0; return 0;	✗
	Overwrite direct flow
	k: a=x; x=0; a=0; return 0;	✓
	Leaking through indirect flow
	l : { if (x) a=0 else a =1}; x=0; return 0;	✗
	Figure 2: Transformed programs
	program p. Definition 2.2 (Security Ordering). A program p leaks more in-
	formation than a program p ′ (written p ≼ p ′ ) if their respective
	Attacker Knowledge are ordered such that
	∀c, m, m ′ .

ACKNOWLEDGMENTS

This work is partially funded by the grant ANR-14-CE28-0014.


David González

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Elías Cueto
email: ecueto@unizar.es

Francesco Dell'

Isola D González

Thermodynamically consistent data-driven computational mechanics

Keywords: Data-driven computational mechanics, GENERIC, Governing equations

come

Introduction

Scientific experiments in large infrastructures are able to produce petabytes of raw data per day. With the irruption of the Internet of Things and the industry 4.0 paradigm, modern factories are also producing breathtaking amounts of data which should be stored, curated and processed so as to identify tendencies and analyze hidden correlations, patterns and, eventually, make decisions in real time. Therefore, being able to predict equations governing physical phenomena from experimental data has utmost importance. A growing interest has arose in the last years in this direction [START_REF] Bongard | Automated reverse engineering of nonlinear dynamical systems[END_REF][START_REF] Schmidt | Distilling free-form natural laws from experimental data[END_REF][START_REF] Battista | Numerical investigation of a particle system compared with first and second gradient continua: deformation and fracture phenomena[END_REF][START_REF] Della Corte | Referential description of the evolution of a 2d swarm of robots interacting with the closer neighbors: perspectives of continuum modeling via higher gradient continua[END_REF][START_REF] Della Corte | Modeling Deformable Bodies Using Discrete Systems with Centroid-Based Propagating Interaction: Fracture and Crack Evolution[END_REF].

There are, however, physical laws of different epistemic nature. One cannot attribute the same importance to the second law of thermodynamics than to the response model of a new material, for instance. Often, these last equations are referred to as constitutive equations and have a phenomenological character. Consequently, the process of distilling physical laws from data in an automated way involves the task of finding relevant invariant structures among them-thus satisfying fundamental conservation equations-but also to be able to identify constitutive equations for the different components of the system under scrutiny.

Under non-equilibrium situations, the problem becomes even harder. An algorithm able to identify the dynamic structure of the system should be able to identify not only the conserved magnitudes (i.e., to find out its Hamiltonian structure), but also its dissipative nature. Many recent attempts in this sense begin by assuming that the (possibly nonlinear) dynamics of the system is of the form

żt = d dt z(t) = f (z(t)), (1)
where z refers to the pertinent set of state variables, needed for an adequate description of the system. Most existing techniques are designed to search among a set of predefined operators (polynomial terms, sinus, cosinus, etc.) those that compone the operator that best fits the experimental results [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF]. Other techniques, such as [START_REF] Daniels | Automated adaptive inference of phenomenological dynamical models[END_REF], progressively add terms, somehow resembling a Taylor series, until an adequate degree of accuracy is reached. More recently, the possibility of operating without any constitutive equation has also been explored [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF][START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF]. In this last work, a set of experimental results substitute constitutive equations by developing a method that operates by finding the experimental result closest to guaranteeing equilibrium. None of the aforementioned techniques, however, impose a priori any fundamental law whose satisfaction is mandatory. Thus, fundamental conservation laws will be satisfied depending on the degree of accuracy of the resulting approximation found by the different algorithms, the noise level in the data, etc. Furthermore, the set of relevant variables z-whose precise meaning will be clear hereafter-could be a priori unknown, or maybe non-relevant variables could have been stored during the experimental campaign. Detecting redundant variables is also an important issue.

In this paper, a method is proposed that guarantees the satisfaction of basic principles. Notably, conservation of energy and positive production of entropy are enforced by construction. This is possible by resorting to the so-called GENERIC structure of the problem [START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF]. The GENERIC paradigm generalizes the Hamiltonian structure of classical mechanics by including the dissipative properties of the problem and seems therefore to be the adequate starting point for this project. In Sect. 2, we review the basics of this paradigm. Then, in Sect. 3 we introduce a strategy for the identification of the GENERIC constituent blocks from experimental data. In Sect. 4, we propose three different examples of the performance of the method that reveals the accuracy of the proposed methodology. Finally, the paper is completed with the discussion of the results.

The GENERIC formalism

Under the framework of non-equilibrium thermodynamics, the GENERIC ("General Equation for Non-Equilibrium Reversible-Irreversible Coupling") formalism establishes a completely general equation of the dynamics of a system under reversible and irreversible conditions, as dictated by the evolution of energy and entropy, respectively [START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF][START_REF] Grmela | Dynamics and thermodynamics of complex fluids. i. development of a general formalism[END_REF]. GENERIC allows for a proper description of the thermodynamical structure of problems within solid and fluid mechanics, thermokinetic effects, chemical reactions or relativistic physics, among other fields in which it has been successfully applied [START_REF] Öttinger | Nonequilibrium thermodynamics: a powerful tool for scientists and engineers[END_REF]. Models constructed upon the GENERIC formalism are guaranteed to posses a sound thermodynamic structure (they preserve the symmetries of the system and therefore guarantee the subsequent conservation/dissipation laws of energy and entropy, respectively, by Noether's theorem), and therefore, it seems to be an appealing choice to start with. In addition, GENERIC paves the way to the development of numerical integration algorithms that can be seen as discrete GENERIC instances and therefore share these interesting properties [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF][START_REF] Romero | Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics: Part i: monolithic integrators and their application to finite strain thermoelasticity[END_REF].

In general, the type of systems we are considering in this work is those that can be described by the laws of classical mechanics. Under this framework, let us consider a closed, isolated system whose evolution is to be described during the time interval I = (0, T ]. The state space of the system will be described by judiciously chosen variables in a space S. The term "judiciously chosen" will be clear hereafter. Therefore, the system will be described at time t ∈ I by a function z t = z(t)

: I → S, z ∈ C 1 (0, T ].
The GENERIC structure of the evolution equations for such a system takes the form

żt = L(z t )∇ E(z t ) + M(z t )∇ S(z t ), z(0) = z 0 . ( 2
)
L is the so-called Poisson matrix and will be responsible for the reversible (Hamiltonian) part of the evolution of the system. E represents the energy of the system, as a function of its particular state at time t, z t . In turn, M represents the friction matrix, responsible for the irreversible part of the evolution of the system. Finally, S represents the entropy of the system for the particular choice of variables z. It is important to note that the actual choice of z variables is irrelevant in the sense that any other complete set of variables (possibly including a different number of them) z = f (z) will produce a valid GENERIC structure of the problem, although probably much more complex. Equation ( 2) is supplemented with the complementary degeneracy conditions, i.e.,

L(z) • ∇S(z) = 0, (3a) M(z) • ∇ E(z) = 0. ( 3b
)
By choosing L skew symmetric and M symmetric, positive semi-definite, one ensures that

Ė(z) = ∇ E(z) • ż = ∇ E(z) • L(z)∇ E(z) + ∇ E(z) • M(z)∇ S(z) = 0, (4)
i.e., one ensures the conservation of energy in closed systems. Equivalently,

Ṡ(z) = ∇s(z) • ż = ∇S(z) • L(z)∇ E(z) + ∇S(z) • M(z)∇ S(z) ≥ 0, (5)
which ensures the fulfillment of the second principle of thermodynamics.

The reason for choosing GENERIC as the starting paradigm for this work is therefore obvious: it guarantees two basic principles such as the conservation of energy and the second principle of thermodynamics. It is under the umbrella of GENERIC that we will guarantee that the found laws will respect basic principles of utmost importance (such as the symmetries of the system, see [START_REF] Oettinger | Beyond Equilibrium Thermodynamics[END_REF]) and that, after this identification, the resulting time integration schemes will show the right conservation properties [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF].

However, it is important to note that the choice of the state space of the system, z, is not obvious. First, the chosen variables must be measurable experimentally, of course. Otherwise, the identification of the system will be complicated. Secondly, and perhaps more importantly, not every possible description of a system can be cast into the GENERIC framework. For systems obeying the classical laws of mechanics, the most detailed (yet often very impractical) level of description is that of molecular dynamics. By successively coarse-graining the system, keeping in the description those variables that evolve in a slow manifold, one arrives, in the limit, at a description in which the sole variables that are involved are those that are conserved in the system along time, i.e., one arrives at the thermodynamics description of the system. During this journey, different descriptions are possible, and not all are suitable for a GENERIC structure [START_REF] Español | Statistical Mechanics of Coarse-Graining[END_REF]. In general, however, GENERIC seems to be the right framework to express a set of equations yet to be found by means of fitting of experimental results. It involves the most general variables of the system, energy and entropy, while ensuring the satisfaction of the most general requirements of thermodynamics.

However, the proposed methodology, which will be clear in the following section, does not look for the identification of a constitutive equation. Instead, what we propose is to work in a purely numerical, discrete way. In [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF], we proposed the concept of constitutive manifold, by assuming that the experimental results (typically twelve-dimensional stress-strain couples for solid mechanics or stress-strain rates for fluid mechanics) live in a certain manifold M whose structure is to be identified, possibly with manifold learning techniques. Once identified, the structure of M will allow for an algorithmic numerical approximation of tangent moduli with which to operate in a purely discrete manner (see also [START_REF] Ibañez | Data-driven nonlinear elasticity. Constitutive manifold construction and problem discretization[END_REF]).

The proposed identification algorithm

It is assumed that the available data will provide the necessary information so as to reconstruct the state of the system at discrete time steps. In other words, we assume to know in advance a set of n meas vector-valued measurements of the already mentioned variables of interest Z = {z 0 , z 1 , . . . , z n meas }. By discretizing Eq. ( 2), one arrives at

z n+1 -z n t = L(z n+1 )DE(z n+1 ) + M(z n+1 )DS(z n+1 ), (6)
where, for simplicity, we denoted z n+1 = z t+ t and where L and M are the discrete version of the Poisson and friction operators, respectively. In turn, DE and DS represent the discrete gradients. In general, matrix L is constant over the process, while matrix M frequently varies.

The proposed governing equation identification algorithm will therefore consist of the identification of the discrete values L and M, on one hand, and DE and DS, on the other. In sharp contrast to previous works, such as [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF], we do not pursue the identification of a particular term (such as polynomials, sinus and cosinus) in the law. Instead, we will approximate the mentioned ingredients L, M, DE and DS numerically and construct with them the already introduced constitutive manifold (see [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF][START_REF] Ibañez | Data-driven nonlinear elasticity. Constitutive manifold construction and problem discretization[END_REF][START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF]).

Within the realm of variational data assimilation, we will assume that the experimental measurements in Z are exact [START_REF] Manzoni | A reduced-order strategy for solving inverse bayesian shape identification problems in physiological flows[END_REF][START_REF] Sullivan | Introduction to Uncertainty Quantification[END_REF]. Of course, data will always be corrupted by noise to some extent, and this will oblige us to consider z in a Bayesian framework. This can be done by stating Eq. ( 2) in the framework of a Fokker-Planck equation, which is always possible [START_REF] Soize | The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions[END_REF]. This is currently the focus of our investigations and will be published elsewhere.

It is worth noting that by identifying the particular form of E(z) we are doing no other thing than identifying what we usually call the "constitutive law" in solid mechanics. In the hyperelasticity framework, we usually denote by ψ or W the strain energy density function and obtain some stress measure (usually the second Piola-Kirchhoff tensor) by means of the computation of the gradient of the elastic energy with respect to some strain measure. In a hyperelastic context, we deal with purely conservative systems, and this will be therefore equivalent to obtaining the Poisson matrix and the ∇ E(z) terms.

Therefore, the proposed algorithm will consist in solving the following (possibly constrained) minimization problem within a time interval J ⊆ I:

μ * = {L, M, DE, DS} = arg min μ ||z(μ) -z meas ||, (7)
with z meas ⊆ Z , a subset of the total available experimental results. Since L, M, DE, DS will in general be dependent on z, the choice of the size of z meas deserves also some analysis: global fits will in general lead to excessive smoothing of the identified results, while using too small experimental sets will lead to ill-posed minimization problems. The general form of matrix L is very easy to find, in general-see the examples that follow-while the structure of M is frequently much more intricate. This is why we prefer, in general, to assume L known (if its form can be determined beforehand) and simplify the problem a bit and to re-formulate it as

μ * = {M, DE, DS} = arg min μ ||z(μ) -z meas ||, (8)
so that only the structure of M and the discrete gradients is sought. This avoids to take into consideration all the constraints related to the skew symmetry of L and improves the convergence of the proposed algorithm. By assuming z approximated in a finite element sense, gradient operators can straightforwardly be cast in matrix form as

DE = Az, (9a) DS = Bz, ( 9b
)
where A and B represent the discrete, matrix form of the gradient operators leading to DE and DS, respectively. Once identified, the manifold structure of {L, M, DE, DS} will also serve to construct, following Eq. ( 6), energy-conserving integration schemes with non-decreasing entropy (see [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF]).

At this point, several identification strategies can be envisaged. These are summarized next.

Identification strategies

To solve problem [START_REF] Daniels | Automated adaptive inference of phenomenological dynamical models[END_REF], several strategies could be implemented, depending on the usage of monolithic or staggered approaches, constrained or not.

Monolithic versus staggered approaches

For some systems, it is possible to easily design experiments under conservative (purely Hamiltonian) conditions. For instance, some systems can be isolated so as to obtain isothermal results with whom to identify the Poisson term of the GENERIC form of the model. Subsequent non-equilibrium experiments will help to identify the dissipative structure of the system. Of course, this is not always possible. One of such examples is that of viscoelastic materials. Viscous and elastic (purely conservative) parts of their behavior cannot be decoupled experimentally. In those cases, the whole set of ingredients of the GENERIC formalism must be identified at once.

Total versus partial identification

We have already mentioned that although Eq. ( 8) incorporates the four ingredients of GENERIC to be identified (i.e., L, M, DE, DS), some of these matrices are easily identified without resorting to any experimental result. For instance, if the system is of mechanical nature and purely Hamiltonian, by choosing z = {q(t), p(t)} (position and momentum of each degree of freedom), the Poisson matrix takes always the form

L(z) = 0 1 -1 0 . ( 10
)
This is the case very often. In other situations, very simple considerations shed light on the form of L and-less frequently-of M.

Constrained or unconstrained minimization

As mentioned before, one of the salient features of GENERIC is the guaranteed conservation of energy and positive entropy production, see Eqs. ( 4), [START_REF] Della Corte | Modeling Deformable Bodies Using Discrete Systems with Centroid-Based Propagating Interaction: Fracture and Crack Evolution[END_REF]. In order to identify the form of the discrete gradient operators, DE, DS, two possibilities arise: to verify a posteriori that the numerically identified terms fulfill these requirements, or to impose it a priori, by solving a constrained minimization procedure:

μ * = {M, DE, DS} = arg min μ ||z(μ) -z meas ||, (11a)
subject to:

L • Bz = 0, (11b) M • Az = 0, (11c)
where restrictions (11b) and (11c) will help, in general, to facilitate the identification procedure while ensuring the thermodynamic consistency of the resulting model.

Results

To validate the proposed methodology, we begin by considering two different finite-dimensional problems whose GENERIC structure is known, so as to be able to determine its degree of accuracy. This allows us to determine the best strategies in order to minimize errors in the identification procedure. Then, an infinitedimensional problem is analyzed that also proves the validity of the suggested strategy.

Simple pendulum

The first problem corresponds to a simple pendulum whose purely Hamiltonian structure, in the absence of any dissipation, responds to Eq. [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF]. Therefore, only the discrete matrix form of the gradient, A, was identified. The accuracy of the proposed identification method was tested by reproducing the dynamics of the system, integrated by constructing the time integrators proposed in [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF], namely the simple scheme outlined by Eq. ( 6). As can be noticed, the agreement is almost perfect. Knowing in advance the form of L, the identification of A employing all the available data (241 time increments) leads to an error vanishing to machine precision (see Fig. 1).

On the contrary, if the identification procedure is applied by taking subsets of five experimental measurements, an error in L2-norm of 0.049% was obtained (see Fig. 2).

Thermoelastic pendulum

As an example of a non-equilibrium problem to start with, we have chosen a problem presented in [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF]. The system consists of two pendula connected by thermoelastic springs. Two masses m 1 and m 2 are thus connected by springs of internal energy e 1 and e 2 and are allowed to oscillate around a fixed point (see Fig. 3). We employ Fig. 2 Results for the pendulum problem. Data (circles) versus numerical integration of the identified system when the identification procedure is applied every five time increments the classical notation of q i , p i , i = 1, 2 for position and momenta, respectively. In turn, springs have entropies s j and longitudes at rest λ 0 j , j = a, b. Under this perspective, the state space of the system can be defined as

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 2 = q 1 }.
The GENERIC structure of the problem will naturally involve the internal energy of the system, as already discussed. It can be computed as the sum of the kinetic energy (coming from the movement of the masses) and the potential energy in the springs, i.e.,

E(z) = K 1 (z) + K 2 (z) + e a (λ a , s a ) + e b (λ b , s b ),
where

λ a = √ q 1 • q 1 , λ b = (q 2 -q 1 ) • (q 2 -q 1 ). m 1 m 2 C a C b Fig. 3 Double thermal pendulum
It is assumed that the temperature in the springs, θ j , changes due to the Joule effect, giving rise to

θ j = ∂e j ∂s j , j = a, b.
Assuming a conductivity κ in the springs, the resulting equations of the system will be

qi = p i m i , ṗi = - ∂ ∂q i (e a + e b ),
ṡ j = κ θ k θ j -1 , with i = 1, 2, j = a, b, k = j.
This gives rise to the following gradients, which are to be identified,

∇ E(z) = f a n a -f b n b , f b n b , p 1 m 1 , p 2 m 2 , θ a , θ b ,
∇ S(z) = (0, 0, 0, 0, 1, 1), with f j , n j , j = a, b, the forces in the springs and their respective unit vector along their direction. Note that we assume them unknown and employ these expressions as ground truth just to test the accuracy of the proposed methodology. Something which is particularly appealing from the structure of GENERIC, already commented in Introduction, is that different choices of state-space variables are equally valid. For instance, one can imagine

S = {z | z = (q 1 , q 2 , p 1 , p 2 , e 1 , e 2 )}, or even S = {z | z = (q 1 , q 2 , p 1 , p 2 , θ 1 , θ 2 )},
and it will produce equivalent results (although with different complexities in the resulting expressions, see [START_REF] Romero | Thermodynamically consistent time-stepping algorithms for non-linear thermomechanical systems[END_REF]).

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , M(z) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 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 κ θ b θ a -κ 0 0 0 0 -κ κ θ a θ b ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
, although, as explained before, we assume M unknown.

The different identification procedures suggested in Sect. 3 allowed to determine the precise form of the gradients of energy and entropy, respectively. Errors in the numerical reconstruction of the pendula trajectories with the just identified terms are reported in Table 1. In that table, the mentioned number of points refers to the size of data sets employed for the identification, i.e., the size of the set z meas in Eq. ( 8).

An example of the reconstructed trajectory of the masses is shown in Fig. 4.

Oldroyd-B fluid

The field of non-Newtonian fluid mechanics is where GENERIC originated and one of the most active in terms of identification of complex constitutive equations. The Oldroyd-B fluid model is a clear example of an equation that can be derived by starting from very different assumptions: from continuum considerations to coarse-graining microscopical descriptions or, finally, by employing a Fokker-Planck probabilistic description [START_REF] Owens | Computational Rheology[END_REF].

The Oldroyd-B model originated as a means to model the flow of diluted polymeric solutions. These polymer chains can be considered both from a purely macroscopic point of view and from a microscopic one, by modeling them as linear dumbbells in a Newtonian medium. From a purely continuum point of view, the constitutive equation of an Oldroyd-B fluid arises from considering the deviatoric part T of the stress tensor σ (the so-called extra-stress tensor) of the form,

T + λ 1 ∇ T = η 0 γ + λ 2 ∇ γ , ( 12
)
where the triangle denotes the nonlinear upper-convected derivative introduced by Oldroyd [START_REF] Walters | The distinctive CFD challenges of computational rheology[END_REF]. The coefficients η 0 , λ 1 and λ 2 are material constants. Finally, following standard notation, γ = (∇ s v) = D represents the strain rate tensor.

If we now consider separately the stress in the solvent (denoted by a subscript s) and polymer (denoted by a subscript p) components as

T = η s γ + τ , this gives rise to τ + λ 1 ∇ τ = η p γ ,
which is the constitutive equation for the elastic stress.

An alternative derivation can be obtained by starting from a population of linear dumbbells immersed in a Newtonian fluid [START_REF] Owens | Computational Rheology[END_REF]. What is here intended is to take the Oldroyd-B model as a paradigm of a class of viscoelastic (and therefore non-Hamiltonian) fluids that encompasses both macroscopic and mesoscopic models. Therefore, we assume that we have a set of pseudo-experimental measurements arising from an Oldroyd-B fluid and try first to identify its constitutive equation in a GENERIC framework and, second, to reconstruct the results by integrating the just found constitutive equations by means of Eq. ( 6).

GENERIC form of the problem

In a GENERIC framework, rather than dealing with constitutive equations for the (elastic) extra-stress tensor, as one intuitively usually starts by, it is often suggested (see, for instance, [START_REF] Pasquali | Consistent scaling of thermal fluctuations in smoothed dissipative particle dynamics[END_REF] and references therein) to start by using the conformation tensor c = r r , as the second moment of the dumbbell end-to-end distance distribution function. In fact, this choice allows for a natural introduction of thermal fluctuations in the model [26]. It can be demonstrated that for modeling viscoelastic fluids, the only variables needed to describe the state of the system are the momentum density p = ρv and the conformation density, ĉ = ρc p c, where v represents the fluid's velocity, ρ is the density and c p is the number of polymer molecules per unit mass (see [START_REF] Mavrantzas | Atomistic monte carlo simulations of polymer melt elasticity: their nonequilibrium thermodynamics generic formulation in a generalized canonical ensemble[END_REF] and references therein). This choice can result somewhat unnatural for scientists used to continuum models. And indeed it is, if one takes into account that the conformation tensor is not easily observable from an experimental point of view. However, the reader can easily notice its similitude to the left Cauchy-Green or Finger strain tensor b = F • F T . Indeed, in a cross-linked rubber, c = N 2 b/3, where N is the number of statistical segments in a polymer chain and its length [START_REF] Pasquali | Consistent scaling of thermal fluctuations in smoothed dissipative particle dynamics[END_REF].

It is well known that the extra-stress tensor takes the form

τ = - n 2ζ 12 ∇ c,
i.e., it is in fact proportional to the convective (or Oldroyd) derivative of the conformation tensor. This magnitude is more suitable as a characteristic internal variable of the fluid than c = r r itself.

Therefore, in order to model a viscoelastic fluid in the GENERIC framework, we must consider, for each of the M fluid particles that compose the (discretized) fluid, their position r j , velocity v j , energy E j and extra-stress tensor τ j . The internal energy E j includes, of course, the kinetic energy contribution from the solvent and the beads, plus the potential energy of solvent-solvent, solvent-bead and bead-bead interaction. Each one of the M fluid particles is assumed to contain j = 1, 2, . . . , N dumbbells and, of course, an entropy S j (E j , τ j , V j ), with V j the volume associated with each particle.

Therefore, using the previous notation, the set of variables employed to describe the phase space of the fluid will be ) for the start-up of a Couette flow (shear flow between two plates, typical of a rheometer). In this case, velocity is imposed at the bottom plate, with h = 1.0 and v shear = 1.0. Notice the overshoot of the velocity profile just before the permanent regime (linear velocity profile)

S = {z = (r j , v j , E j , τ j , j = 1, 2, . . . , M) ∈ (R 2 × R 2 × R × R 4 ) M }. (13)

Implementation details

Although the just introduced strategy is entirely equivalent to the existing data-driven simulation approaches by Ortiz and coworkers [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF][START_REF] Kirchdoerfer | Data driven computing with noisy material data sets[END_REF], in the sense that it performs a constrained minimization problem at every integration point, alternative approaches exist. Instead of fitting point-wise the best possible constitutive values provided by data, it is possible to do it in a non-local manner, by employing model order reduction approaches. Model reduction consists, essentially, in finding the slow manifold in which the system evolves. Therefore, for a (full order) problem formulated under the GENERIC paradigm, ensuring that its reduced-order version also verifies the GENERIC structure is not obvious. The interested reader can consult [START_REF] Öttinger | Preservation of thermodynamic structure in model reduction[END_REF]. The obvious advantage is the drastic reduction in the number of degrees of freedom.

In this case, we employed a very simple strategy based on the employ of the proper orthogonal decomposition (POD) of the pseudo-experimental results [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Probability theory[END_REF][START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF]. An alternative, perhaps more accurate, is to employ some form of nonlinear dimensionality reduction, or manifold learning strategy, as in [START_REF] Millán | Nonlinear manifold learning for model reduction in finite elastodynamics[END_REF].

Our pseudo-experimental results were generated with a CONNFFESSIT approach [START_REF] Laso | Calculation of viscoelastic flow using molecular models: the CONNFFESSIT approach[END_REF][START_REF] Cueto | Meshless stochastic simulation of micro macro kinetic theory models[END_REF] to the problem of a Couette flow of an Oldroyd-B fluid. The Couette flow is the pure shear flow between two rigid plates produced typically in rheometers, and thus the most frequent form of experimental results for this type of fluids. The CONNFFESSIT approach combines finite elements and stochastic techniques for the solution of micro-macro-flows of viscoelastic flows.

Given the rheometric character of the flow (in other words, its homogeneous character), the code discretizes only the one-dimensional distance h between the rheometer plates, by employing one hundred finite elements (and, consequently, 101 nodes). We consider the start-up of a flow with a Reynolds number of 0.1 and Weissenberg number 0.5. Each finite element node is equipped with 10,000 dumbbells.

A simple proper orthogonal decomposition (POD) of the 150 snapshots of the velocity profile (see Fig. 5) revealed that the flow can be expressed in terms of a few degrees of freedom (see Fig. 6). In fact, only six modes produce a decay on the eigenvalues of more than four orders of magnitude. Therefore, instead of the 101 nodes employed to build the mesh, our identification procedure employed the just mentioned six degrees of freedom.

Thus, for each of the six retained degrees of freedom, the Poisson and friction matrices are known to have the form

L k (z) = ⎛ ⎜ ⎝ 0 1 0 0 -1 0 1 -1 0 -1 0 0 0 1 0 0 ⎞ ⎟ ⎠ , M k (z) = ⎛ ⎜ ⎝ 0 0 0 0 0 1 1 0 0 1 1 0 0 0 0 1 ⎞ ⎟ ⎠ , k = 1, . . . , n POD ,
with n POD the number of employed POD modes. They are, however, assumed unknown, and their form is determined by the proposed method. In this case, the identification includes both L and M without any significant difficulty. Note that, due to the particular form of the Couette flow, the variables in S (see Eq. ( 13)) turn out to be one-dimensional and hence the 4 × 4 size of the Poisson and friction matrices.

With these ingredients, the pseudo-experimental and identified results for nodes 1, 11, . . . , 101 along the vertical axis of the flow are represented in Fig. 7.

As in previous examples, the accuracy of the proposed methods depends on the number of time samples employed to perform the minimization problem. In this case, by performing a constrained minimization (thus enforcing conservation of energy and strictly positive entropy dissipation) with sets of 150 snapshots, the error in the velocity field turned out to be of 10 -8 . The unconstrained minimization strategy, on the contrary, was unable to find any suitable solution for the problem.

Numerical integration of an arbitrary flow

Of course, the final goal of this procedure is not to reproduce the experimental results but to be able to integrate different flow regimes once the constitutive manifold of the particular fluid has been identified. Assume, for simplicity, that we have identified the discrete gradient operators at time t n , DE = A n z n and DS = B n z n , and the Poisson and friction matrices L and M, respectively-that turned out to be constant for this problem-for two flows under different conditions, say Weissenberg, W e, numbers. We have, therefore, the values L, M and the triplets A 1 n , B 1 n , z 1 n and A 2 n , B 2 n , z 2 n , together with the current state of the system, z n . The objective is then to numerically integrate, without any explicit knowledge on the constitutive equations of an Oldroyd-B fluid, the state equations for the flow under a third W e value, characterized by A obj , B obj . Assume known the value of the weights w 1 and w 2 that will make Eq. ( 6) hold under the new flow conditions, namely,

z obj n+1 = z n + t • L A obj + M B obj z n = z n + w 1 t • L A 1 n + M B 1 n z 1 n +z n + w 2 t • L A 2 n + M B 2 n z 2 n
with L and M also known. Some simple algebra leads to the following conditions:

A obj z n = w 1 A 1 n z 1 n + w 2 A 2 n z 2 n , B obj z n = w 1 B 1 n z 1 n + w 2 B 2 n z 2 n .
Computing the interpolating weights w 1 , w 2 , . . . , w , with the number of neighboring experimental results to interpolate with, can be done, for instance, by resorting to nonlinear dimensionality reduction techniques such as locally linear embedding (LLE), which produces these weights as a by-product of the dimensionality reduction process [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF].

By doing so, the results of integrating a Couette flow at W e = 0.75 from the experimental results obtained at W e = 0.5 and W e = 1.0 (thus, w 1 = w 2 = 0.5 in a straightforward manner, for this simple illustration) are shown in Fig. 8. In it, the data-driven results (we recall again that no constitutive equation was employed to perform the integration) are compared toward the result of a direct numerical integration. Results for the stress component τ xy are shown. L2-norm error of the velocity field shown to be 0.4988% in this particular case.

Poiseuille flow

Obviously, the aim of the just developed technique is not to be able to reproduce experimental results. Rather, it is to be able to simulate situations for which no experimental result exist. Usually, rheologists employ rheometers to fully characterize the constitutive equations of a particular fluid. The simplest of these rheometers provokes simple shear (Couette) flows, as in the previous sections. In this section, we will try to demonstrate that, once the constitutive manifold has been properly characterize by means of the necessary experimental results, the approach just developed could be employed for the simulation of flows for which no previous result exists.

So, let us assume that we have performed enough experiments so as to characterize the constitutive manifold. In this case, we have simulated a big amount of shear (Couette) flows to characterize the fluid. With the help of these (pseudo-)experimental results, we are about to simulate the start-up of a Poiseuille flow. By taking a time increment of t = 0.006, we interpolate from the set of available (pseudo-)experimental results the A and B matrices that are employed to advance the flow in time. Figure 9 represents the evolution of the τ xy component along time, for both the reference solution and the just computed one. For this flow, the mean relative error for this value is between 4.13 and 16.14% for the different nodes of the discretization (which is actually one-dimensional). Of course, these errors can be lowered by employing finer meshes and smaller time steps.

Within the realm of data-intensive computational mechanics, this paper presents a novel methodology able to integrate the equations of motion without resorting to any constitutive equation (the so-called data-driven computational mechanics framework first established in [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF]). The main novelty in the approach herein consist in the development of a thermodynamically consistent approach to the problem that guarantees the fulfillment of the first and second principles of thermodynamics, independently of the quality of the experimental results. This is in sharp contrast to existing approaches in the literature, either the ones working without any constitutive assumption (see [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF][START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF][START_REF] Ibañez | Data-driven nonlinear elasticity. Constitutive manifold construction and problem discretization[END_REF][START_REF] Lopez | A manifold learning approach for integrated computational materials engineering[END_REF]) or those that identify the form of the governing equations by choosing among a predefined set of equations, [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Peherstorfer | Data-driven operator inference for nonintrusive projection-based model reduction[END_REF].

The present approach is based on the usage of the GENERIC paradigm, a framework that guarantees the strict energy conservation and positive entropy production, along with the symmetries of the system, and that generalizes the Hamiltonian form of classical mechanics to non-equilibrium systems. The proposed methodology thus identifies the proper GENERIC structure of the problem. The method identifies either the whole GENERIC structure composed by the Poisson and friction matrices of the problem, plus the precise numerical form of the gradients of the energy and entropy. However, for some complex problems it has been found that the form of the Hamiltonian and dissipative structures of the problem is easy to identify and that the method works considerably better by just finding the form of the gradients of energy and entropy, respectively. This is equivalent, of course, to finding in a numerical, point-wise manner, the form of the constitutive manifold that the authors defined in a previous work (see [START_REF] Ibañez | A manifold learning approach to data-driven computational elasticity and inelasticity[END_REF]).

But the proposed methodology not only provides with a suitable method to unveil the precise form of the Hamiltonian and dissipative structures of the problem. It also provides a robust and thermodynamically consistent way to integrate these equations without employing any constitutive equation (the so-called datadriven computational mechanics paradigm).

The aforementioned properties have been demonstrated by solving different problems: from simple pendulums with discrete degrees of freedom to complex viscoelastic fluids. In this last case, the proposed technique has been combined with a classical model order reduction technique (POD) that allows to highly simplify the problem and to greatly diminish the number of degrees of freedom and variables to identify.

Fig. 1

 1 Fig.1Results for the pendulum problem. Data (circles) versus numerical integration of the identified system

Fig. 4

 4 Fig.[START_REF] Della Corte | Referential description of the evolution of a 2d swarm of robots interacting with the closer neighbors: perspectives of continuum modeling via higher gradient continua[END_REF] Results for the thermal pendulum problem. Trajectories of the two masses. Data (circles) versus numerical integration of the identified system. Monolithic identification on sets of |z meas | = 34 data points

Fig. 5

 5 Fig.[START_REF] Della Corte | Modeling Deformable Bodies Using Discrete Systems with Centroid-Based Propagating Interaction: Fracture and Crack Evolution[END_REF] Velocity profile at every t = 0.01 seconds (t ∈ (0, 1.5]) for the start-up of a Couette flow (shear flow between two plates, typical of a rheometer). In this case, velocity is imposed at the bottom plate, with h = 1.0 and v shear = 1.0. Notice the overshoot of the velocity profile just before the permanent regime (linear velocity profile)

Fig. 6

 6 Fig. 6 POD modes and related eigenvalues for the Couette flow

Fig. 7 Fig. 8 Fig. 9

 789 Fig. 7 Experimental versus reconstructed phase space variables for the Oldroyd-B fluid. a Displacement, b velocity, c energy and d x y component of the extra-stress tensor, τ xy . Every variable is represented for nodes 1, 11, . . . , 101 along the height of the channel, h

Table 1

 1 Errors (L2-norm) for the double thermal pendulum

			Strategy
			Staggered			Monolithic
			84 points (%)	17 points (%)	167 points (%)	34 points (%)
	Constrained		0.000001	0.0684		0.0023		0.3072
	Unconstrained		0.000018	0.3114		0.00036		0.0399
		14
						Data-ball 1
		12				Simulation-ball 1 Data-ball 2
						Simulation-ball 2
		10
		8
	Y	6
		4
		2
		0
		-2 -2	0	2	4	6	8	10
					X

Acknowledgements This work has been supported by the Spanish Ministry of Economy and Competitiveness through Grants number DPI2017-85139-C2-1-R and DPI2015-72365-EXP and by the Regional Government of Aragon and the European Social Fund, research group T88.