A minimizing-movements approach to GENERIC systems

We present a new time discretization scheme adapted to the structure of GENERIC systems. The scheme is variational in nature and is based on a conditional incremental minimization. The GENERIC structure of the scheme provides stability and convergence of the scheme. We prove that the scheme can be rigorously implemented in the case of the damped harmonic oscillator. Numerical evidence is collected, illustrating the performance of the method.


Introduction
The aim of this note is to discuss a new variational time-discretization scheme adapted to the structure of General Equations for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC).Introduced by Grmela & Öttinger, this formulation provides a unified frame for describing the time evolution of physical systems out of equilibrium in presence of reversible and irreversible dynamics [36].
Let y denote the state of a closed, nonequilibrium physical system and let E(y) and S(y) be the corresponding total energy and total entropy, respectively.The GENERIC formulation of the time evolution of the system reads y ′ = L(y) DE(y) + K(y) DS(y). (1.1) Here, y ′ denotes the time derivative, L is the antisymmetric Poisson operator, and K is the symmetric and positive definite Onsager operator (details in Section 2).The gist of the GENERIC system (1.1) is that conservative and dissipative dynamics are clearly separated.Under a specific compatibility condition, see (2.4) below, this entails that a solution t → y(t) to (1.1) conserves energy and accumulates entropy in a quantifiable way, namely d dt E(y) = 0 and d dt S(y) = DS(y), K(y)DS ≥ 0.
The conservation of energy and the quantified dissipative character are the distinguishing traits of the GENERIC system (1.1).To replicate these properties at the discrete level leads to so-called structure-preserving approximations.These have drawn interest in the last years, giving rise to different numerical solutions adapted to different applicative contexts.
Numerical schemes conserving energy can be found for instance in [16,40,[42][43][44]47]; see [8] for a review and [49] for a contribution explicitly focusing on GENERIC formulations.These schemes are often discrete-gradient methods, where gradients of functionals are specifically modified in order to fulfill a discrete chain rule and exactly replicate conservation [16,17,19,30].In the thermomechanical context, structure-preserving discretizations either in terms of the absolute temperature [9,10,15,38] or of the internal energy or the entropy [5,6] have been obtained.The reader is referred to [28] for an approach to open systems.
GENERIC integrators able to conserve energy and accurately describe entropy accumulation have been proposed in [37], where however some limitations are also mentioned.In particular, energy and Onsager operators have to be suitably modified in a time-step dependent manner in order energy conservation to hold.Integrators are actually constructed in case of single dissipation mechanism only (were K be a matrix, it would have to have rank one) and no convergence theory is provided.The discretization of the damped harmonic oscillator is addressed in [37], where nonetheless the temperature is given by a prescribed heat bath.In this case, explicit solutions have to be used in order to specify the GENERIC integrator.
To the best of our knowledge, all available numerical schemes directly target GENERIC systems in their differential form (1.1).Our focus is here on a timediscretization of variational nature instead, fitting into the general scheme of socalled minimizing movements [1,2].In particular, at all discrete steps we aim at solving a specific minimization problem.We start by defining the entropyproduction potential Ψ(y, ξ) = ξ, K(y)ξ /2 and denoting by Ψ * (y, •) its conjugate in the second variable.Then, the GENERIC relation (1.1) in [0, T ] can be equivalently rewritten in terms of the scalar equation for all t ∈ [0, T ].The reformulation of dissipative systems in terms of scalar equations as (1.2) is usually referred to as De Giorgi's Energy-Dissipation Principle [12,29,33].This principle has already been applied in a variety of different contexts, including generalized gradient flows [4,46], curves of maximal slope in metric spaces [2,11], rate-independent systems [34,41], and optimal control [39].For GENERIC systems, the reformulation of (1.1) in terms of such variational principle has been already presented in [14].
Our new minimizing-movements approach to (1.1) consists in tackling a discrete version of relation (1.2).Assume to be given a time partition 0 = t 0 < t 1 < • • • < t N = T with steps τ i = t i − t i−1 , and an initial state y 0 .We define the time-discrete trajectory {y i } N i=0 by letting y 0 = y 0 and subsequently perform the minimization min y −S(y) for i = 1, . . ., N .The latter is nothing but a localized and discretized version of the De Giorgi's Energy-Dissipation Principle (1.2).The minimizing-movements scheme above has been introduced in [25].The theory in [25] is however restricted to the case of state-independent operators L and K, which severely limits the applicability to real GENERIC systems.In addition, the convergence analysis in [25] relies on a suitable set of a priori assumptions, leaving open the discussion whether these can be met in practice.
The aim of this note is then threefold.First, we extend the reach of the numerical method to include the case of state dependent operators L(y) and K(y), hence covering the full extent of the GENERIC theory (Section 2).Our main result is the conditional convergence of Theorem 2.1.Second, we provide the detailed analysis of a specific case, the damped harmonic oscillator, in which the above-mentioned convergence assumptions can actually be proved to hold (Section 3).Eventually, we present numerical experiments assessing the performance of the minimizingmovements scheme and compare it to the classical implicit Euler scheme (Subsection 3.5).

The minimizing-movements scheme for GENERIC systems
In this section, we recall the structure of a GENERIC system [36] by specifying the assumptions on functionals and operators that will be used throughout.Moreover, we formulate our minimizing-movements scheme and present a conditional convergence result, namely Theorem 2.1.
The GENERIC system y ′ = L(y) DE(y) + K(y) ∂S(y) a.e. in [0, T ] (2.1) is defined by specifying the quintuple (Y, E, S, L, K).In the following, the state space Y is assumed to be a reflexive Banach space.The functionals E and S represent the total energy and the total entropy of the system, respectively.We assume E to be Fréchet differentiable, with strongly-weakly continuous differential DE, and −S : Y → (−∞, ∞] to be proper and lower semicontinuous with singlevalued and strongly-weakly continuous Fréchet subdifferential ∂(−S).Recall that one has ξ ∈ ∂(−S)(y) iff S(y) > −∞ and In the following, we also make use of the obvious notation −∂S = ∂(−S).
The operators L and K define a Poisson and an Onsager geometric structure on Y , respectively.In particular, for all states y ∈ Y we assume that L(y) and K(y) are linear and continuous from Y * to Y .Moreover, L(y) is required to be antisymmetric, L * (y) = −L(y), and to fulfill the Jacobi identity {{g 1 , g 2 }, g 3 } + {{g 2 , g 3 }, g 1 } + {{g 3 , g 1 }, g 2 } = 0. Here, g i denotes any differentiable function on Y and the Poisson bracket is defined as {g, g}(y) = Dg(y), L(y)Dg(y) , where •, • denotes the duality pairing between Y * (dual) and Y .Moreover, we assume the strong-weak continuity (2.2) The mapping K(y) is asked to be symmetric and positive definite, namely K(y) = K * (y) ≥ 0. We associate to K the so-called entropy-production potential Ψ : Y × Y * → [0, ∞) given by Ψ(y, ξ) = ξ, K(y)ξ /2 and let Ψ * be its conjugate in the second variable, namely Ψ * (y, η) = sup ξ ( ξ, η − Ψ(y, ξ)).We also assume the lower semicontinuity of the sum of the entropy-production potential and its dual, that is (2.3) In addition, functionals and operators are asked to satisfy the crucial noninteraction condition (2.4) This condition ensures that the system of dissipative actions K(y) ∂S(y) does not spoil energy conservation and the system of conservative actions L(y) DE(y) does not contribute to dissipation.Indeed, by assuming sufficient smoothness one checks that = ∂S(y), K(y)∂S(y) ≥ 0. (2.5) The noninteraction condition (2.4) hence implies that trajectories y solving (2.1) have constant total energy and entropy rate ∂S(y), K(y)∂S(y) .In particular, the entropy is nondecreasing and entropy production results solely from irreversible processes.In computing (2.5) we have used the chain rule (d/dt)S(y) = ∂S(y), y ′ almost everywhere.This is classical in case −S is (a regular perturbation of) a convex functional [7,Prop. 3.3,p. 73].The reader is referred to [45] for a general discussion out of the convex case.Before moving on, let us remark that the structure of GENERIC is geometric in nature.Indeed, it is invariant by coordinate changes.Let Then, the quintuple ( Ỹ , Ẽ, S, L, K) satisfies the above structural assumptions and the GENERIC structure (2.1) can be rewritten as ỹ′ = L(ỹ) D Ẽ(ỹ) + K(ỹ) ∂ S(ỹ).We now reconsider the discussion leading to (1.2) and specify it further by remarking that relation (2.1) is actually equivalent to the inequality − S(y(t)) for all t ∈ [0, T ].The equivalence between (2.1) and (2.6) follows from Fenchel's relations.Applied to the entropy-production potential Ψ, these relations read where the subdifferential is taken in the second variable only.By noting that ∂Ψ(y, ∂S(y)) = K(y) ∂S(y) and using (2.7)-(2.8),one can prove the equivalences (2.8) (2.4) (2.7) In particular, the last left-to-right implication follows by integration while the rightto-left counterpart from the nonnegativity of the integrand, given (2.4) and (2.7).The minimizing-movements scheme corresponds to a discretization of inequality (2.6).To each time partition 0 = t 0 < t 1 < • • • < t N = T , we associate the time steps τ i = t i − t i−1 and the diameter τ = max τ i .Given the vector {y i } N i=0 ∈ Y N +1 , we introduce the backward piecewise constant and piecewise linear interpolations We define the incremental functional by By letting y 0 = y 0 we find the discrete solution {y i } N i=0 by subsequently solving the minimization problem min y G(τ i , y i−1 ; y) for i = 1, . . ., N. (2.9) Note that, for all τ > 0 and y i−1 ∈ Y with S(y i−1 ) > −∞, the map y → G(τ i , y i−1 ; y) is strongly lower semicontinuous because of the lower semicontinuity of −S, the lower semicontinuity (2.3) of Ψ * + Ψ, the weak-strong continuity of DE and ∂S, and the continuity (2.2) of L. In order to solve problem (2.9) one has hence to check that y → G(τ i , y i−1 ; y) is strongly coercive.
The main result of this section is the following conditional convergence result.
Theorem 2.1 (Conditional convergence).Under the above assumptions, let a sequence of partitions 0 = t n 0 < • • • < t n N n = T be given with τ n → 0 as n → ∞, and let {y n i } N n i=1 be such that y n 0 = y 0 .Assume that ) ) Then, up to a subsequence, y n ⇀ y in H 1 (0, T ; H), where y is a solution of the GENERIC system (2.1) in the sense of inequality (2.6).
The conditional convergence result of Theorem 2.1 relies on the possibility of solving the inequality G(τ n i , y n i−1 ; y n i ) ≤ 0 up to a small, controllable error, and establishing some a priori bounds on the discrete solution.The validity of these conditions has to be checked on the specific problem at hand.In the coming Section 3 we give an example of a situation where (2.10)-(2.13)actually hold.
In case −S is convex, an example of {y n i } N n i=0 fulfilling (2.10) are the solutions of the implicit Euler scheme whenever available.Indeed, such where the last equality follows from (2.8) and (2.18).Note that the incremental minimization problem (2.9) may have no solutions, even if the Euler scheme is solvable.In the purely dissipative case L = 0, the existence of a solution to (2.9) is ensured as soon as τ * > 0 exists, so that, for all τ ∈ (0, τ * ) and y i−1 ∈ Y with S(y i−1 ) > −∞, the function y → −S(y) + τ Ψ(y, (y − y i−1 )/τ ) is strongly coercive.The latter follows whenever −S has strongly compact sublevels and would imply in particular that the Euler scheme is also solvable.We anticipate however that the example discussed in Section 3 is not purely dissipative and −S, albeit convex, does not have strongly compact sublevels.
Let us mention that, in specific applications, the bounds (2.11)-(2.13)may follow from (2.10).For instance, this would be the case if the coercivity holds for some c > 0 and all y, η ∈ Y and ξ ∈ Y * , namely in case K(y) is positive definite and bounded, uniformly with respect to y.This however does not apply to the example of the damped harmonic oscillator from Section 3, for K is singular there.
The quadratic nature of the entropy-production potential could be generalized to the case of p-growth with p > 1 without any specific intricacy.In particular, one could consider the polynomial case ∂Ψ(y, ξ) = K(y) ξ p−2 Y * ξ and coercivity (2.20) would then read Y * ∂S(y), K(y) DS(y) ≥ 0. Let us further remark that the conditional convergence result of Theorem 2.1 can serve as an a-posteriori tool to check the convergence of time-discrete approximations {y n i }, regardless of the method used to generate them.In particular, relation (2.10) can be seen as a sort of a-posteriori error estimator.
We conclude this section by pointing out that the statement of Theorem 2.1 would hold also for other minimizing-movements schemes, where nonlinearities are possibly evaluated explicitly.In particular, instead of G one could consider the functional G where y j , j = 1, 2, 3, are independently chosen to be either y i or y i−1 (or else, for instance (y i +y i−1 )/2).Under conditions (2.10)-(2.13),convergence would again follow as in Theorem 2.1.

The minimizing-movements scheme for the damped harmonic oscillator
In this section, we analyze the simplest GENERIC system fulfilling conditions (2.10)-(2.13).We consider the case of the damped harmonic oscillator, namely ) Here, q ∈ R represents the position of the harmonic oscillator and θ > 0 is its absolute temperature.Note that the case of θ being the given, constant temperature of a surrounding heat bath is considered in [37] instead.The positive constants m, ν, κ, λ and c are the mass of the oscillator, the viscosity of the medium, the elastic modulus, the thermal-exchange coefficient, and the heat capacity, respectively.By letting p ∈ R be the momentum of the harmonic oscillator, we rewrite (3.1)-(3.2) as the first-order system ) System (3.3)-(3.5)can be written in the GENERIC form (2.1) by letting y = (q, p, θ) ∈ Y = R 2 × (0, ∞) represent the state of the harmonic oscillator and by defining the total energy and the total entropy as and S(q, p, θ) = −λq + c + c ln θ.
In particular, S = −∂ θ Φ and E = p 2 /(2m) + Φ + θS.Note that both E and S are smooth, −S is convex, and S(y) > −∞ iff θ > 0. The operators L, K : R 3 → R 3×3 are given by One readily checks that L is antisymmetric and a tedious computation shows that the Jacobi identity holds.On the other hand, K is symmetric and positive semidefinite, while not being invertible.Note that K(y) has rank two for all p = 0, so that the construction in [37] does not apply here.
We devote the remainder of this section to prove that the convergence result of Theorem 2.1 applies to the scheme (3.6)-(3.9),namely that conditions (2.10)-(2.13)are fulfilled.In particular, we check that the minimization problem admits a solution {y i } N i=0 = {(q i , p i , θ i )} N i=0 (Subsection 3.1), that such solution is unique if time steps are small enough (Subsection 3.2), that each such solutions fulfills condition (2.10) in the stronger sense G(τ i , y i−1 ; y i ) ≤ 0 (Subsection 3.3), and that the bounds (2.11)-(2.13)can be established, independently of the partition (Subsection 3.4).Eventually, numerical simulations are presented in Subsection 3.5.
By using the constraints (3.8)-(3.9),we can express q i and θ i in terms of p i and the values (q i−1 , p i−1 , θ i−1 ) in the form where the function f is defined by By substituting the two expressions in the minimum problem (3.6)-(3.9),we can reduce it to a minimization in the variable p i only.Indeed, problem (3.6)-(3.9) is equivalent to min pi where we have defined Again, as F is smooth on its domain, in order to solve the minimization problem (3.11) we just need to prove that F has bounded sublevels.This follows from the fact that F (p i ) → ∞ if f (p i ) → 0 + and that f (p i ) > 0 iff p i belongs to a bounded interval, depending on (q i−1 , p i−1 , θ i−1 ).In fact, as f is quadratic and Hence, for all given (q i−1 , p i−1 , θ i−1 ) ∈ R 2 × (0, ∞), we can find a solution p i of (3.11).We conclude from (3.10) that (q i , p i , θ i ) solves (3.6)-(3.9).In particular, we have that θ i > 0.

Uniqueness of minimizers.
Let us remark that the minimizers (q i , p i , θ i ) of (3.6)-(3.9)need not be unique in general.Still, the function F is strictly convex for τ i small enough, depending on the values (q i−1 , p i−1 , θ i−1 ) and on material parameters.
Let us start by rewriting F as where h i = (−p i−1 + τ i κq i−1 )/(1 + τ 2 i κ/m) and r i is the affine function Using the following facts and differentiating F in (3.12) twice, we obtain This condition depends on the size of the time-step and on material parameters and implies that F is strictly convex, thus admitting a unique minimizer p i .Consequently, problem (3.6)-(3.9)has a unique minimizer (q i , p i , θ i ).
We can make the time-step dependent bound (3.13) on τ i more explicit by directly computing the maximum of f : In Subsection 3.4 we will check that above right-hand side can be bounded by a positive constant c 0 in terms of the initial data and of material parameters only, uniformly with respect to i, see (3.22).Condition (3.13) can hence be strengthened by asking the partition to be such that This stronger condition ensures the unique solvability of the minimization problem (3.6)-(3.9)for all i = 1, . . ., N .
We may hence rewrite equation (3.17) as with γ i , δ i , and ǫ i depending on α i , β i , θ i−1 , and data.Note that i κ) 2 < 0. We infer from θ i−1 > 0 that ǫ i > 0. Hence, the second-order equation (3.19) has a unique solution θ i > 0. This uniquely defines p i and q i via (3.18) and (3.15), respectively.In particular, the implicit Euler scheme (3.15)-(3.17)has a unique solution y e i = (q i , p i , θ i ).
Let now y i be a solution to the incremental minimization problem (3.6).Owing to (2.19), we conclude that = 0, where we have also used the fact that y e i fulfills the constraints (3.7)-(3.9).In particular, condition (2.10) holds in the even stronger form G(τ i , y i−1 ; y i ) ≤ 0 for i = 1, . . ., N. (3.20)

3.4.
A priori bounds.We now prove that condition (2.11)-(2.13) of Theorem 2.1 hold.This will follow by checking a priori bounds on minimizers {y i } N i=0 of (3.6), independently from the time partition.
Let us start by remarking that the energy is nonincreasing.By (3.8), one can equivalently rewrite the constraint (3.9) as A rearrangement of the terms leads to Hence, E(y i ) is nonincreasing and E(y i ) = E(y i−1 ) iff q i = q i−1 and p i = p i−1 .By taking the sum for i = 1, . . ., j for j ≤ N , we find that where we have set The nonnegative terms D q j and D p j exactly quantify the energy dissipativity of the scheme.Owing to (3.21), we obtain the uniform bound where, here and in the following, the symbol C denotes a generic positive constant depending on the initial data y 0 and on material parameters, but not on the time partition.Bound (3.22) and constraint (3.8) imply that Moreover, we readily check that Let us take the sum for i = 1, . . ., j for j ≤ N in (3.20) getting −S(y j ) + ψ * j + ψ j ≤ −S(y 0 ), where As ψ * j and ψ j are nonnegative, we infer that S(y i ) is nondecreasing.In particular, −c ln θ j ≤ −S(y 0 ) − λq j + c ≤ C ∀j = 1, . . ., N.
and consequently, for some given θ min depending just on y 0 and the material parameters.This ensures that |DS(y It follows from the bound (3.22) on q i and θ i that S(y i ) = −λq i + c + c ln θ i ≤ C for all i = 1, . . ., N .We conclude that ψ N ≤ −S(y 0 ) + S(y N ) ≤ C and Using the fact that θ i is uniformly bounded by (3.22), we can hence bound Eventually, we infer from constraint (3.9) that where (q, p, θ) solves the damped harmonic oscillator system (3.3)-(3.5)with initial value (q(0), p(0), θ(0)) = (q 0 , p 0 , θ 0 ).Note that solutions of (3.3)-(3.5)are unique.Hence, convergence (3.30) holds for the whole sequence of discrete solutions, not just for a subsequence.As E and S are smooth, y n is bounded, and θ i ≥ θ min > 0, we find that E(y(•)) → E(y(•)) ≡ E(y 0 ) and S(y(•)) → S(y(•)) uniformly, together with their time derivatives of all orders.More precisely, upon remarking that ≤ Cτ n , we may use (3.21) and check that In particular, we control the energy dissipation as follows: 3.5.Numerical tests.We record here some numerical evidence on the performance of the minimizing-movements scheme (3.6)-(3.9).All computations are performed in Matlab.In the following, we choose Figures 1-3 illustrate the numerical reference solution and the time-discrete solutions for τ = 1/4.Both the minimizing-movements and the Euler scheme dissipate energy, see Figure 3 left.This is of course an undesired effect, which is however attenuated as τ → 0, see (3.31).Energy dissipation seems to be more pronounced for the Euler scheme.On the other hand, entropy is nondecreasing for the minimizingmovements scheme whereas the Euler scheme shows a nonmonotone entropy, which is nonphysical, see Figure 3  we compute the uniform errors of the temperature component of the solution and of the energy with respect to the numerical reference solution.
As τ converges to 0, our computations confirm that both the minimizing-movements and the Euler scheme are of order τ , see also (3.31).Let us mention that a proof of first-order convergence for the minimizing-movements scheme in the nondissipative regime (L = 0, K independent of the state) is given in [25,Prop. 4.3].Both schemes do not conserve energy.Still, the minimizing-movements scheme is more accurate than the Euler scheme when the time steps are large.Position q (left) and momentum p (right) with respect to time for the numerical reference solution (red, dotted), the minimizing-movements scheme (solid), and the Euler scheme (blue, dash-dotted).

Figure 1 .
Figure1.Position q (left) and momentum p (right) with respect to time for the numerical reference solution (red, dotted), the minimizing-movements scheme (solid), and the Euler scheme (blue, dash-dotted).

Figure 2 .
Figure2.Temperature with respect to time for the numerical reference solution (red, dotted), the minimizing-movements scheme (solid), and the Euler scheme (blue, dash-dotted).

Figure 3 .
Figure 3. Energy (left) and entropy (right) as function of time for the numerical reference solution (red dotted), the minimizingmovements scheme (solid), and the Euler scheme (blue, dashdotted).