Structure-preserving integrators for dissipative systems based on reversible-irreversible splitting

We study the optimal design of numerical integrators for dissipative systems, for which there exists an underlying thermodynamic structure known as GENERIC (general equation for the nonequilibrium reversible-irreversible coupling). We present a frame-work to construct structure-preserving integrators by splitting the system into reversible and irreversible dynamics. The reversible part, which is often degenerate and reduces to a Hamiltonian form on its symplectic leaves, is solved by using a symplectic method (e.g., Verlet) with degenerate variables being left unchanged, for which an associated modified Hamiltonian (and subsequently a modified energy) in the form of a series expansion can be obtained by using backward error analysis. The modified energy is then used to construct a modified friction matrix associated with the irreversible part in such a way that a modified degeneracy condition is satisfied. The modified irreversible dynamics can be further solved by an explicit midpoint method if not exactly solvable. Our findings are verified by various numerical experiments, demonstrating the superiority of structure-preserving integrators over alternative schemes in terms of not only the accuracy control of both energy conservation and entropy production but also the preservation of the conformal symplectic structure in the case of linearly damped systems.


Introduction
As an introduction to this article on structure-preserving integrators for dissipative systems, we first summarize the state-of-the-art of the literature and then provide a description of the GENERIC formulation and its properties. The introduction ends with an outline of the article.

State-of-the-art in structure-preserving integrators
In the last few decades, considerable effort has been devoted to developing structure-preserving integrators for Hamiltonian systems. It has been demonstrated that the so-called symplectic integrators, which preserve the symplectic structure, have superior long time behavior compared to their nonsymplectic counterparts, and should be preferred in practice [16,30,32]. On the other hand, there has been growing interest in designing appropriate numerical methods for gradient flows [2,15,22,51,62] that respect their underlying properties. In contrast to the symplectic structure, the conformal symplectic structure [5,6,10,20,38,46] for Hamiltonian systems that are perturbed by a linear damping (which can be thought of as a special case of the Rayleigh dissipation) has been less studied. It is also worth mentioning that variational integrators [23] and specialized Runge-Kutta methods [21] have also been used to solve dissipative systems. It turns out that thermodynamically admissible evolution equations for nonequilibrium systems have a more general (including an additional variable known as entropy) and well-defined structure known as GENERIC [13,47,49,50], which possesses the following distinct features: (i) conservation of the total energy; (ii) separation of the reversible and irreversible dynamics; (iii) the reversible dynamics preserves a Poisson structure; (iv) entropy production is unaffected by the reversible dynamics; (v) nonnegative entropy production rate.

GENERIC formulation
The GENERIC formulation of the time evolution for nonequilibrium systems is given by where x is the set of independent variables required to describe a given nonequilibrium system, E and S represent, respectively, the total energy and entropy as functions of the independent variables x, and L and M denote the antisymmetric Poisson matrix and the positive semidefinite (symmetric) friction matrix, respectively. Note that both L and M can also depend on the independent variables x so that the fundamental time evolution equation (1) could be highly nonlinear. We also point out that ∂/∂x in (1) simply implies the partial derivative although it typically denotes the functional derivative when x is a function/field. Moreover, (1) is supplemented by two degeneracy conditions: and Equations (2)-(3) indicate the conservation of the entropy by the reversible dynamics (i.e., the L contribution) and the conservation of the total energy in a closed system by the irreversible dynamics (i.e., the M contribution), respectively. Note that "reversible" and "irreversible" dynamics (in thermodynamics) are simply the names of the two fundamental contributions to the time evolution equation (1), and should not be confused with similar terms in other subjects. The rank of M has the interpretation of the number of dissipative processes taking place in the system. (See more discussions on the formulation of the GENERIC framework in [13,47,49,50].) The usefulness and maturity of the GENERIC framework have been illustrated in a very large number of successful applications in a wide range of areas in Appendix E of [47] (see also a most recent review of [48] and references therein). In particular, despite its simple form, we believe that the irreversible dynamics in (1) is the most general form of meaningful irreversible equations in nonequilibrium thermodynamics-it is a belief based on both a very large variety of successful examples and statistical mechanics, so that it can be called knowledge (in particular, as this belief is widely accepted in the nonequilibrium thermodynamics community).
In order to further demonstrate the general properties of L and M , the respective Poisson and dissipative brackets are often adopted: where A and B are sufficiently regular (and real-valued) functions of the independent variables x. With the help of the two brackets and the chain rule, the time evolution equation of an arbitrary function A can then be written as More specifically, the Poisson bracket (4) inherits the antisymmetry of L, and satisfies the Leibniz rule, {AB, C} = A{B, C} + B{A, C} , where C is another arbitrary sufficiently regular (and real-valued) function of the independent variables x. In addition, the Poisson bracket is required to satisfy the Jacobi identity, The dissipative bracket (5) inherits the symmetry of M , and also satisfies the Leibniz rule, The positive semidefinite nature of M leads to the nonnegativeness condition which implies the second law of nonequilibrium thermodynamics (i.e., the entropy production rate is always nonnegative), This article addresses the long-standing challenge of how to preserve the underlying structures when numerically discretizing GENERIC systems in practice. Although in recent years, this topic has attracted increasing attention [26,27,44,52], to the best of our knowledge, there are no such numerical integrators in the literature. Unlike common approaches that are based on exact energy conservation, we propose in this article a framework to construct structurepreserving integrators for dissipative systems, i.e., GENERIC integrators (also known as metriplectic integrators [12,24,25,42,43] in the mathematical literature), based on splitting the reversible and irreversible dynamics. The topic of structure-preserving integrators for GENERIC/metriplectic systems is the counterpart and generalization of the theory of symplectic integrators for Hamiltonian systems.

Outline of the article
The rest of the article is organized as follows. We give specific definitions of GENERIC integrators and discuss their requirements in numerical discretizations in Section 2. In Section 3, we propose a framework to construct split GENERIC integrators based on reversible and irreversible splitting, the generality of the framework is demonstrated in examples of linearly damped systems in Section 3.1 as well as in a more challenging (and fully coupled) case of two gas containers exchanging heat and volume in Section 3.2. Section 4 presents various numerical experiments to investigate the performance of the two split GENERIC integrators introduced in this article. Our findings are summarized in Section 5.

Definitions of GENERIC integrators
In this section, we provide the definitions of GENERIC integrators and discuss their requirements when numerically discretizing a system in practice.

Full GENERIC integrators
We recall the definition of (full) GENERIC integrators given in [49]. Analogous to the definition of symplectic integrators for Hamiltonian dynamics [45], a mapping, x 0 → x h , is said to be a full GENERIC integrator if it corresponds to a continuous time evolution of a modified whereẼ h andM h represent the modified energy and friction matrix associated with the integrator, respectively, satisfying a modified degeneracy condition: That is, given initial conditions x(0) = x 0 , the analytical solution of (14), x(t), should agree with what we obtain from the integrator at time h, i.e., x(h) = x h . A full GENERIC integrator x → x h , which can be thought of as the formal solution of (14), possesses the following structure: Similar to symplectic integrators for Hamiltonian dynamics, the modified energy,Ẽ h , is strictly conserved by a GENERIC integrator. The physical energy E is expected to remain close to the modified energy,Ẽ h , even for long integration periods. Additionally, the modified friction matrix,M h , should not introduce any additional dissipative processes not present in the original matrix M . We point out that full GENERIC integrators may only be available in special cases, for instance, a full GENERIC integrator in the case of a damped harmonic oscillator, where analytical solutions of the GENERIC system can be obtained, was proposed and discussed in [49]. However, it should be noted that it is highly unlikely that analytical solutions would be available for general GENERIC systems. (Nevertheless, it might be eventually possible to recognise a full GENERIC integrator without exact solutions.) Therefore, in what follows we introduce a framework to construct "split" GENERIC integrators.

Split GENERIC integrators
Inspired by recent developments on splitting methods [1,[28][29][30][31][33][34][35]59], we consider to split the reversible and irreversible parts of the GENERIC system in such a way that the reversible dynamics, which is often degenerate but possesses a Hamiltonian form on its symplectic leaves, can be integrated by using a symplectic method (e.g., Verlet) with degenerate variables being left unchanged, while the irreversible part (gradient flow) can be solved in such a way that as many structure elements as possible can be preserved (see more references on the challenging task of structure preservation on manifolds in [2,15,22,36,37,51,62]). An interesting question for the split GENERIC integrators is: under what conditions do a modified energy and an associated friction matrix, satisfying the modified degeneracy condition (15), exist? If they exist, how much do we know about their respective forms? GENERIC integrators share some common features of GENERIC systems discussed at the beginning of this article, which can also be thought of as the requirements for GENERIC integrators. Denoting the Jacobian matrix of the independent variables x as Ω, we have (i) preservation of the Poisson structure for the reversible dynamics: (ii) nonnegative entropy production rate: S(x h ) ≥ S(x 0 ); (iii) the modified degeneracy condition (15) is satisfied with the other (2) being unchanged; (iv) preservation of the rank of the friction matrix: rank(M h ) = rank(M ).
Note that the satisfaction of the modified degeneracy condition (15) may be based on a truncated modified energy as discussed in Section 3.1.2. As pointed out in [53], it has been proved in [67] that there cannot exist an integrator for "non-integrable" Hamiltonian dynamics that preserves both the symplectic (Poisson) structure and the energy (Hamiltonian). In fact, it has been discussed in [61] that the preservation of either property has its advantages and disadvantages. While previous attempts to construct structure-preserving integrators for dissipative systems have been relying on the exact conservation of energy (i.e., the energyconserving discrete gradient methods [8,40,54], see more discussions in Section 4.1.2), there is no obvious reason why integrators that preserve the Poisson structure for the reversible dynamics should be ignored.

Construction of split GENERIC integrators based on reversibleirreversible splitting
In this section, we discuss the construction of GENERIC integrators based on splitting the reversible and irreversible parts of the system. In order to satisfy the modified degeneracy condition (15), we explore the possibility of adjusting the irreversible part using a modified friction matrix that corresponds to a modified energy associated with the symplectic integrator used for the reversible part.

Linearly damped systems
We first consider a linearly damped system that possesses a natural GENERIC structure (1) with independent variables x = (q, p, S), where q and p represent the position and momentum of the particle, respectively, and S is the entropy of the surrounding thermal bath. While S is an independent variable and thus ∂S/∂x = (0, 0, 1), the total energy of the GENERIC system is given by where H(q, p) represents the Hamiltonian of the particle, U (q) denotes the potential energy, and T S is the energy of the thermal bath. Given the antisymmetric Poisson matrix and the positive semidefinite (symmetric) friction matrix where constant parameters m, γ, and T represent the mass of the particle, the damping rate, and the constant temperature of the thermal bath, respectively, the equations of motion of the GENERIC system can be written aṡ where F (q) = −U ′ (q) is the conservative force. Note that in this particular case the symplectic leaves are given by the (q, p) subsystem within the reversible dynamics for constant entropy S.

The YBABY method
Following the discussions in Section 2.2, we suggest to split the GENERIC system (20)- (22) into reversible and irreversible parts, Moreover, we can always use a symplectic method (e.g., Verlet) for the reversible dynamics on its symplectic leaves (this is possible in the setting of linearly damped systems (20)- (22) where S is an independent variable), while for linearly damped systems with the total energy (17) the irreversible dynamics is exactly solvable (with q being left unchanged) Therefore, we can apply the Verlet method to integrate the reversible part (24) and then further split the exact solver (26)- (27), e hL Y , for the irreversible part (25) to composite a symmetric splitting method, termed "YBABY", as where exp (hL f ) denotes the phase space propagator associated with the corresponding vector field f , with L f being the corresponding generator. The generators for each part of the GENERIC system may be written out as follows: Thus, the generator for the GENERIC system can be written as The integration steps of the YBABY method read: The order of convergence of a splitting method can be determined by using the Baker-Campbell-Hausdorff formula [16,30,32]. For general operators A and B, we have where with A, B = AB − BA being the commutator. Subsequently, we can work out where Therefore, a symmetric splitting typically gives second order convergence whereas a nonsymmetric one is generally first order. One can then obtain the associated operator of the YBABY method:L which indicates formally second order convergence for the YBABY method (29). Note that the order of convergence can also be demonstrated by using the Taylor series expansion for the solutions, but the procedure is often tedious. Note also that in principle higher order methods can also be constructed, as in Hamiltonian dynamics [66], by suitably composing the operators.
We would also like to point out that while all three subsystems can be solved exactly in linearly damped systems, in cases where the irreversible part is not exactly solvable (see the example of two gas containers exchanging heat and volume in Section 3.2) it is important to solve the irreversible part by using a numerical method that is at least second order so that an overall second order convergence is expected. Alternatively, one could solve the irreversible part by using a numerical method, which could be first order (e.g., the Euler method), and its adjoint method for half a step each, it can be shown that the resulting YBABY † method is self-adjoint (or symmetric) and typically has even order (see more discussions in [16,30,32]). However, such a method could become implicit, for instance the adjoint method of the Euler method is the implicit backward Euler method.
In the case of γ = 0, the YBABY method reduces to the Verlet method with degenerate variable S being constant, which is a well-known symplectic method that preserves the Poisson structure for the reversible dynamics [16,30,32]. Therefore, in order to guarantee the preservation of the Poisson structure for the reversible dynamics, in what follows we will apply the Verlet method for the reversible part, unless otherwise stated.
For linearly damped systems, it has been demonstrated in [5,6] that numerical methods that preserve the underlying "conformal symplectic" structure [38] are advantageous over alternative schemes. Moreover, high order conformal symplectic and ergodic schemes for stochastic Langevin equation have also been investigated [20]. Definition 3.1. A numerical method is said to be conformal symplectic if the symplectic two form decays exponentially with a constant decay rate, i.e., where ∧ represents the wedge product and K > 0 is the constant decay rate. Similarly, a numerical method is said to be symplectic if the symplectic two form is preserved, i.e., dq h ∧ dp h = dq ∧ dp .
We point out that if the prefactor in front of dq ∧ dp is initially not in an exponential form, we can always rewrite it into an exponential form as long as it is a constant value between zero and one. Following [5,20], we can show that the YBABY method (29) is conformal symplectic: dq n+1 ∧ dp n+1 = e −γh/2 dq n+1 ∧ dp n+2/4 , = e −γh/2 dq n ∧ dp n+2/4 , = e −γh dq n ∧ dp n .
in which case the decay rate is the physical damping rate, i.e., K = γ.
We have so far verified the second order convergence for the YBABY method, and its preservation of the Poisson structure for the reversible dynamics as well as the conformal symplecticity. However, it is unclear under what conditions there exist a modified energy and an associated friction matrix as in (15). To this end, in what follows we modify the irreversible part of the system as discussed at the beginning of this section.

The mYBABY method
It is well known that if a symplectic method is used for the reversible dynamics (24), there exists a modified Hamiltonian,H h , in the form of a (typically infinite) series expansion obtained by using backward error analysis [55], which is exactly preserved by the symplectic integrator [16,30,32]. In the example of the Verlet method, the modified Hamiltonian is given byH In order to identify a modified energy conserved by a GENERIC integrator, we can replace the original energy E (17) by a modified energy,Ẽ h =H h + T S, and then try to explore whether we can construct an associated friction matrix,M h , in such a way that the modified degeneracy condition (15) is satisfied. However, it is unlikely that we can find such a friction matrix due to the infiniteness of the series expansion (and often complicated higher order terms) in the modified energy. Nevertheless, we can truncate the series expansion of the modified energy to certain order in practice, which will introduce some perturbations to the modified energy. For instance, we can use the Verlet method for the reversible part, and then truncate the modified energy up to second order, introducing a perturbation of order four to the modified energy, to obtaiñ Subsequently, we can construct the associated modified friction matrix in the fashion of backward error analysis [16,32,55]: whereỹ h is assumed to be a truncated series expansion up to second order with In order to satisfy the modified degeneracy conditioñ which leads toỹ the following condition has to be satisfied which has a solution Thus, the modified friction matrix can be written as where the "modifying factor" is given by Moreover, the modified friction matrix induces a small (second order) perturbation of the physical entropy production As a result, the irreversible part, incorporating the modified friction matrix (56), becomes which can be solved exactly (with q being left unchanged) In this case, the generator for the modified irreversible dynamics becomes By replacing the Y piece by Y m in the YBABY method (29), we can similarly define a symmetric splitting method, termed "Y m BABY m " or "mYBABY", as where the associated operator can be worked out by applying the Baker-Campbell-Hausdorff formula [16,30,32] asL which indicates formally second order convergence for the mYBABY method (63). It can be easily shown that all four requirements listed in Section 2.2 are satisfied for the mYBABY method. The integration steps of the mYBABY method read: Note that in the case of the "modifying factor" (57) being unity, the mYBABY method (63) reduces exactly to the YBABY method (29). It can be shown that the truncated energyẼ h (49) is the truncated modified energy, up to second order, for the mYBABY method (63), based on the fact that: (i) the Verlet method for the reversible dynamics preservesẼ h (49) at second order; (ii) the exact solver for the irreversible dynamics preservesẼ h (49) exactly. In principle, we could truncate the modified energyẼ h at higher orders (e.g., fourth, sixth,. . . ) than that of (49), which would lead to higher orders for the overall methods if the irreversible dynamics can be solved exactly. Moreover, it might be more appropriate to refer those GENERIC integrators that incorporate the truncation of the modified energy to "pseudo-GENERIC integrators" (in a sense similar to pseudo-symplectic integrators that preserve the symplectic structure only to certain orders [3]).
It can be further shown that the mYBABY method (63) preserves the conformal symplectic structure if the Hessian of the potential energy is a constant, i.e., U ′′ (q) = C. That is, following (47), we have dq n+1 ∧ dp n+1 = e −γmh/2 dq n+1 ∧ dp n+2/4 , = e −γmh/2 dq n ∧ dp n+2/4 , = e −γmh dq n ∧ dp n , where which can be thought of as a modified decay rate compared to the damping rate in the YBABY method (47). Note that the preservation of the conformal symplectic structure is, in the literature, often associated with a decay rate of exactly the damping rate as in the YBABY method. Therefore, we may interpret that the YBABY method preserves the conformal symplectic structure in a "strong" sense whereas the mYBABY method preserves the conformal symplectic structure in a "weak" sense.

Two gas containers exchanging heat and volume
In order to demonstrate the generality of our framework introduced in Section 3.1.2, we also consider an example of two (ideal) gas containers exchanging heat and volume (see Fig. 1 and Exercises 3 & 9 in [47] for more details) with independent variables x = (q, p, S 1 , S 2 ), where q and p, respectively, represent the position and momentum of the separating wall of mass m, while S 1 and S 2 are, respectively, the entropies of the two subsystems. In this case the total energy is given by where E 1 and E 2 are, respectively, the internal energies of the two subsystems with the following relationships to their associated entropies and volumes (i.e., the Sackur-Tetrode equation for ideal gases) where k B is the Boltzmann constant,ĉ is another constant that is needed to ensure the argument of the logarithm dimensionless, and it is assumed that the two subsystems contain the same number of particles, N . The volumes of the two subsystems are given by where A c is the area of the cross section and 2L g is the length of the container. Given the antisymmetric Poisson matrix and the positive semidefinite (symmetric) friction matrix where the positive constant parameter α determines the strength of the heat exchange, and T 1 and T 2 are, respectively, the temperatures of the two subsystems, related to the associated internal energies by the resulting equations of motion of the GENERIC system can be written aṡ Although the motion of the wall is assumed to be frictionless (i.e., there is no explicit damping term as in linearly damped systems in Section 3.1), all oscillations of the separating wall have to be damped since they induce (time dependent) temperature differences and thus a heat flux with entropy production. Alternatively, an analysis of the equations (80)-(83) linearized around equilibria indicates that the system would always relax to equilibrium. We would also like to point out that, unlike linearly damped systems considered in Section 3.1 where the (q, p) dynamics may be viewed as being independent of the entropy, the (q, p) dynamics in this case "strongly" depends on the dynamics of (S 1 , S 2 ), and vice versa-it is a fully coupled GENERIC system. As in Section 3.1.1, we could also split the system (80)-(83) into reversible and irreversible parts: for which we can use the Verlet method for the reversible dynamics, with degenerate variables S 1 and S 2 being constants, while a suitable method can be used to solve the irreversible dynamics. We would again like to construct a modified energy and an associated friction matrix as in (15). To this end, following the procedures in Section 3.1.2 we first identify the modified energy associated with the Verlet method used for the reversible dynamics based on which we can subsequently work out the derivatives of the truncated modified energy up to second order: and We can then construct the associated modified friction matrix in the fashion of backward error analysis as where the modifying factors are given by and respectively. The modified irreversible part, incorporating the modified friction matrix (88), is now given by The modified friction matrix again induces a small (second order) perturbation of the physical entropy production Both YBABY and mYBABY methods are similarly defined in this setting as for linearly damped systems in Section 3.1. However, we are unable to solve the modified irreversible dynamics exactly here, thus a second order explicit midpoint method is suggested to approximate the modified irreversible dynamics (91) while the Verlet method is still used for the reversible dynamics, with degenerate variables S 1 and S 2 being constants. Overall, the two split GENERIC integrators are both expected to be second order.
We would like to point out that in some cases it might be beneficial to replace the explicit midpoint method by alternative (higher order and/or higher accuracy) methods. Moreover, inspired by the subsampling techniques popular in large-scale Bayesian sampling [34,60], it might be computationally highly advantageous (especially in high dimension) to decompose the positive semidefinite modified friction matrix into nonoverlapping principal submatrices (a principal submatrix can be obtained by selecting a subset of rows and the same subset of columns) that are still positive semidefinite. Having avoided directly solving a high dimensional gradient flow, we could instead solve each of the decomposed and much smaller subsystems with a significantly reduced computational overhead (even with high accuracy). A thorough investigation of this direction is beyond the scope of this article, and will be left for future work.
It is also worth mentioning that when the modified irreversible dynamics has to be approximated by certain numerical methods, the truncated modified energy is expected to be preserved in an "approximation" sense. A detailed analysis of the effect of the approximation is also beyond the scope of this article, and will be left for future work.

Numerical experiments
In this section, we conduct various numerical experiments to examine the performance of the two split GENERIC integrators introduced in this article.

Simulation details
In the case of linearly damped systems, we consider one-dimensional examples of a damped harmonic oscillator (i.e., U (q) = kq 2 /2), for which an analytical solution can be obtained [49], as well as a damped nonlinear oscillator (i.e., U (q) = −k cos(q)) where the argument of the cosine function should be dimensionless and this is achieved by fixing the unit of length via the initial position q 0 . The equations of motion of both linearly damped systems can be simplified by dimensional analysis [64]. Without loss of generality, in both cases we choose the basic units (mass, time, temperature, and length, respectively) as m = k = T = 1 and q 0 = 2, where the initial position was particularly chosen to demonstrate the nonlinear effects in the damped nonlinear oscillator. Subsequently, the equations of motion of both linearly damped systems involve only the single dimensionless parameter of γ ≥ 0. Moreover, in both cases we chose p 0 = 0 as more general values of the initial momentum essentially correspond to a shift of the initial time. Since we are more interested in the deviation from the initial entropy than its absolute value, we set the initial entropy to be zero in both cases.
In the other case of two gas containers (where A c = L 2 g ) exchanging heat and volume, we chose the basic units of mass and length as m = L g = 1, respectively. We further set N k B = 1, which fixes a characteristic macroscopic unit of entropy, the counterpart of T = 1 (i.e, the thermodynamic unit) in the previous examples of linearly damped systems. In order to fix the fourth unit of time, α = 0.5 was chosen so that: (i) the period of the oscillation is of order one; (ii) there are enough oscillations to collect statistical data (larger values of α lead to faster decay of the amplitude of the oscillation). Furthermore, initial conditions of (q, p, E 1 , E 2 ) = (1, 2, 2, 2) were used (i.e., the separating wall is initially in the middle of the container with an initial velocity).
In all three cases, the positions appeared to be oscillating with the amplitudes decaying exponentially. The total simulation time T s in each case was thus chosen so that t = T s is the time at which the amplitude of the oscillation was reduced to approximately 1/e times its initial value.
Denoting h as the integration stepsize and subsequentlyN = T s /h the number of integration steps, the root-mean-square error (RMSE) of observable φ is defined as follows: whereφ i and φ i represent the numerical approximation at time ih and its corresponding exact (reference) value, respectively. In order to demonstrate the superiority of structure-preserving integrators over alternative schemes, we compare the two split GENERIC integrators introduced in this article with the explicit third order Runge-Kutta (RK3) method used also in [5] as well as the average discrete gradient (ADG) method [17]. The choice of the RK3 method is clearly arbitrary, while other methods are typically second order, it serves as a good example of a higher order method that is not structure-preserving.

The third order Runge-Kutta method
Rewriting GENERIC systems in a compact form asẋ(t) = f (t, x) with initial conditions x(0) = x 0 , the RK3 method is given by where with t n = nh, n = 0, 1, 2, . . . . Note that the RK3 method is neither symplectic nor conformal symplectic. Thus, it does not preserve the Poisson structure for the reversible dynamics. It is also worth mentioning that the two split GENERIC integrators introduced in this article at each step typically require only one force calculation, which often dominates the computational cost per step especially for large-scale simulations, whereas three force calculations are needed for the RK3 method.

The average discrete gradient method
The so-called discrete gradient methods [8,40,54], which are also known as discrete derivative methods [11], have often been used for the time integration of dissipative systems [14,19,26,27,52,[56][57][58]63]. For instance, they have recently been suggested to temporally discretize the Landau collision operator in an attempt to preserve its metriplectic/GENERIC structure [25]. However, as stated in [7,39], discrete gradient methods are generally not symplectic for the symplectic leaves and thus the Poisson structure of the reversible dynamics is not preserved. Therefore, those discrete gradient methods do not belong to either of the GENERIC integrators defined in Section 2. Moreover, discrete gradient methods are typically implicit, in which case iterative methods (e.g., Newton's method) are needed to approximate the solutions at each step. Therefore, discrete gradient methods could be considerably more time-consuming than alternative explicit methods depending on not only the stopping criterion for the iterating procedure [7] but also the size of the linear system that needs to be solved at each iteration. However, in the special case of a damped harmonic oscillator (i.e., U (q) = kq 2 /2), we can work out the integration steps without the iterating procedure. More precisely, we rewrite the GENERIC system (20)- (22) as which is discretized by where the matrixS(x n , x n+1 ) approaches S(x) in the limits of x n+1 → x n and h → 0, while the discrete gradient∇E(x n , x n+1 ) satisfies the following conditions: We consider a midpoint discretization forS, i.e., and the ADG [17] for∇E, i.e., in which case the ADG method reduces to the implicit midpoint method, which is second order and symplectic (for the symplectic leaf with γ = 0) [16,18] q n+1 = q n + hp n+1/2 /m , where q n+1/2 = q n + q n+1 /2 and p n+1/2 = p n + p n+1 /2. One might be surprised how, with the irreversible dynamics being switched off (i.e., γ = 0), the energy-conserving ADG method (or the implicit midpoint method) can also be symplectic for the symplectic leaves, which seems to "contradict" the findings of [67] (see discussions in Section 2.2). However, we point out that in the case of a harmonic oscillator, the corresponding Hamiltonian subsystem is in fact integrable, in such a special case the ADG method preserves not only the energy but also the Poisson structure. Moreover, we can easily solve (104)-(105) to obtain and subsequently (if 4m + h 2 k > 2mhγ) dq n+1 ∧ dp n+1 = e −γ ADG h dq n ∧ dp n , where We can see from (109)-(110) that the "symplectic two form" of the ADG method decays exponentially with a constant decay rate, thus the ADG method in this special case preserves the conformal symplectic structure in a "weak" sense. Furthermore, when the irreversible dynamics is switched off (i.e., γ = 0), (109) reduces to dq n+1 ∧ dp n+1 = dq n ∧ dp n , which indicates that the ADG method preserves the Poisson structure for the reversible dynamics. We emphasize here that if the force is nonlinear or alternative discrete gradient approximations (e.g., the midpoint discrete gradient [11], which was used in several methods compared in [26]) are used, the preservation of the conformal symplectic structure and the Poisson structure for the reversible dynamics is expected to be violated while the iterating procedure seems to be unavoidable, which could result in a substantial computational overhead. The ADG method is unsurprisingly implicit in both cases of the damped nonlinear oscillator (i.e., U (q) = −k cos(q)) and two gas containers exchanging heat and volume. While the former is similar to the damped harmonic oscillator case except replacing kq in (98) by k sin(q), whose ADG is still analytically integrable, the latter is more involved. To be more precise, we rewrite the GENERIC system (80)-(83) as where with the constantB being defined asB In this case, the ADG for∇E q (x n , x n+1 ) is no longer analytically integrable, and thus approximated by using the trapezoidal rulē As a result of the approximation, the exact total energy conservation of the ADG method is expected to be violated. (Note that one may rewrite the GENERIC system (80)-(83) with independent variables x = (q, p, E 1 , E 2 ). However, the same issue of violating the exact total energy conservation for the ADG method is still expected.)

Damped harmonic oscillator
We first consider the damped harmonic oscillator (i.e., U (q) = kq 2 /2) example, where the analytical solution is available using the same set of parameters (except q 0 ) in [49]. It is of great importance that numerical approximations of GENERIC systems have (i) a good conservation of the total energy and (ii) a faithful production of the physical entropy, both of which were compared in Fig. 2. We compare the performance of the two split GENERIC integrators  (17) (left) and entropy (right) against stepsize by comparing the two split GENERIC integrators introduced in this article with the third order Runge-Kutta (RK3) method and the average discrete gradient (ADG) method, which conserves the total energy exactly (i.e., up to machine precision) and thus is only included in comparisons of the entropy production, with a damping rate of γ = 0.01 and a total simulation time of T s = 200 in a standard setting of a damped harmonic oscillator as described in Section 4.1. The stepsizes tested began at h = 0.0094 and were increased incrementally by 30% until around h = 0.5. Note that, with this set of parameters, the damped harmonic oscillator is "underdamped" (i.e., the position of the particle oscillates around zero with the amplitude exponentially decreasing to zero) and the associated period is T p ≈ 2π. Dashed black lines represent the second and third order convergence as indicated.
with that of the RK3 method and the ADG method. Since the ADG method conserves the total energy exactly (i.e., up to machine precision) in this setting [53], its results will not be included in comparisons of the energy conservation. According to the dashed order lines, the RK3 method shows third order convergence whereas other methods are all second order as expected. Among the second order methods, the mYBABY method in all the cases we tested outperforms the YBABY method in terms of the RMSE in both quantities. Particularly, in terms of the entropy production, mYBABY is remarkably one order of magnitude more accurate than YBABY, which is slightly outperformed by the ADG method. In both cases, despite its higher computational overhead, the higher order RK3 method is only more accurate than either of the split GENERIC integrators when the stepsize is relatively small, especially for the mYBABY method. The evolutions of the absolute error in the total energy from various methods were compared (against the exact value of E 0 = 2) and plotted in Fig. 3. We can see from the figure that, with a stepsize of h = 0.1 (left panel), the absolute error of the RK3 method rises quickly before eventually settling down, while the absolute errors of the two split GENERIC integrators oscillate strongly with the amplitudes decreasing. Consistent with our findings in Fig. 2, the absolute error of mYBABY is largely smaller than that of YBABY. Nevertheless, both mYBABY and YBABY methods are more accurate than the RK3 method with a relatively large stepsize. The behavior is rather similar with a larger stepsize of h = 0.5 (right panel) except the magnitude of the error obtained from each method is considerably larger than that with a smaller stepsize. Figure 4 compares the control of the entropy production from a variety of methods. The behavior of the YBABY, mYBABY, and RK3 methods are similar to that of Fig. 3. Interestingly, with a stepsize of h = 0.1 (left panel), the absolute error of the ADG method, while oscillating, initially grows before decreasing while the absolute error of mYBABY, also oscillating, is constantly smaller than that of ADG. The behavior is again very similar with a larger stepsize of h = 0.5 (right panel) except the magnitude of the errors.
We also compare in Fig. 5 the decay of the oscillation amplitude represented by the "local maximum" (in logarithm) of the numerical solution of the position, which characterizes the preservation of the conformal symplectic structure. It can be seen from the figure that while the decay rate of the YBABY method is preserved (almost indistinguishable from the reference decay rate of the damping rate γ), the RK3 method, which is not conformal symplectic, exhibits a clear drift. It can be also observed (and verified) that both mYBABY and ADG decay at slightly different rates of γ m (73) and γ ADG (110), respectively, compared to the reference decay. This indicates that both mYBABY and ADG in this particular case preserve the conformal symplectic structure in a "weak" sense (see discussions at the end of Section 3.1.2).

Damped nonlinear oscillator
We also investigate the performance of various methods with the damped nonlinear oscillator (i.e., U (q) = −k cos(q)), where the reference solution was obtained by using the RK3 method with a very small stepsize of h = 0.001. It turns out that the performance of those methods is very similar to that in the case of the damped harmonic oscillator in Section 4.2. Therefore, we only present the results of the accuracy control of both energy conservation and entropy production as in Fig. 2. According to the dashed order lines in Fig. 6, the RK3 method again exhibits third order convergence whereas the other methods appear to be second order. Moreover, we can see from both panels of the figure that the mYBABY method again comfortably outperforms the YBABY method. The RK3 method, with a higher computational overhead, is again in both cases only more accurate than either of the split GENERIC integrators when the stepsize is relatively small. As mentioned in Section 4.1.2, the time-consuming iterating procedure had to be adopted for the ADG method in this nonlinear case. Since the ADG method conserves the total energy up to machine precision, we only include it for comparisons of the entropy production. As can be seen from the right panel of Fig. 6 that the ADG method is more accurate than the YBABY method but is outperformed by the mYBABY method despite its higher computational overhead.

Two gas containers
We further examine the performance of various methods in the case of two gas containers exchanging heat and volume described in Section 3.2, where the reference solution was again obtained by using the RK3 method with a very small stepsize of h = 0.001. While the RK3 method still shows third order convergence, the other methods appear to be second order as expected, according to the dashed order lines in Fig. 7. The performance of the two split GENERIC integrators and the RK3 method is largely similar to that in the previous two examples. More precisely, in both cases the mYBABY method still clearly outperforms the YBABY method, and the two split GENERIC integrators appear to be more accurate than the RK3 method unless the stepsize is relatively small. Unlike the damped nonlinear oscillator example in Section 4.3, the ADG is not analytically integrable and had to be approximated, resulting in errors in the total energy. As a result, we can see from Fig. 7 that the performance of the ADG method is almost indistinguishable from that of the YBABY method in terms of the accuracy control of the energy conservation on the left panel, while the former is outperformed by the latter in terms of the accuracy control of the entropy production on the right panel. In both cases, the ADG method is clearly outperformed by the mYBABY method, although the latter only preserves the truncated modified energy in an "approximation" sense. Moreover, we would like to point out that while the two split GENERIC integrators are both explicit, the implicit ADG method is computationally much more time-consuming (in this particular case, the evolution of the system was obtained by using the iterative Newton's method at each step in which a linear system associated with a 4 × 4 Jacobian matrix was repeatedly solved).

Conclusions
We have given specific definitions of GENERIC integrators that preserve the underlying thermodynamic structures. In order to construct such integrators, we have presented a framework by splitting a GENERIC system into reversible and irreversible parts. The former, which is often degenerate and reduces to a Hamiltonian form on its symplectic leaves, is solved by a symplectic (Verlet) method (with degenerate variables being left unchanged) for which an associated modified Hamiltonian (and subsequently a modified energy) can be obtained by using backward error analysis. The modified energy is subsequently used to construct a modified friction matrix associated with the irreversible part in such a way that the modified degeneracy condition (15) is satisfied. Following the framework, the mYBABY method has been proposed, which, along with another split GENERIC integrator of the YBABY method, is expected to be second order and typically require only one force calculation at each step. Between the two split GENERIC integrators, we have observed that mYBABY clearly outperforms YBABY in all the cases tested, indicating the importance of satisfying the modified degeneracy condition (15). We have demonstrated by conducting a variety of numerical experiments (including linearly damped systems and two gas containers exchanging heat and volume) that, in terms of the accuracy control of both energy conservation and entropy production, the two split GENERIC integrators (particularly the mYBABY method) are more accurate than the higher order RK3 method unless the stepsizes are relatively small, not to mention the latter requires three force calculations at each step. While the two split GENERIC integrators preserve the conformal symplectic structure for linearly damped systems, RK3 fails and exhibits a clear drift in the decay of the oscillation amplitude of the numerical solutions.
Since the ADG method conserves the total energy up to machine precision, we do not include it in comparisons of the energy conservation for linearly damped systems. It turns out that in both examples of linearly damped systems, the ADG method appears to be more accurate than YBABY, but (despite the use of the time-consuming iterating procedure in the damped nonlinear oscillator case) outperformed by mYBABY in terms of the accuracy control of the entropy production. The ADG is not analytically integrable and had to be approximated in the case of two gas containers exchanging heat and volume, leading to errors in the total energy. As a result of that approximation, the ADG method is as accurate as the YBABY method in terms of the accuracy control of the energy conservation, while the former is outperformed by the latter in terms of the accuracy control of the entropy production. In both cases, although preserving the truncated modified energy in an "approximation" sense, the mYBABY method is clearly more accurate than the ADG method. This indicates that in cases where approximations have to be made in the ADG method, it could lose its "builtin" advantage of exact conservation of the total energy and be outperformed by alternative methods (especially mYBABY). Moreover, we would like to emphasize again that the implicit ADG method is considerably more time-consuming than the two split GENERIC integrators due to the use of the iterative Newton's method where a linear system associated with a 4 × 4 Jacobian matrix was repeatedly solved at each step. It is anticipated that the computational overhead of discrete gradient methods could be substantially increased for large systems, making them unfavorable compared to explicit structure-preserving integrators (especially mYBABY) in practice.
It is worth mentioning that it might be possible to design GENERIC integrators without an explicit construction of the modified energyẼ h . This is related to the question whether the irreversible dynamics (i.e., a vector field) is "compatible" with the canonical transformation associated with the time step h. Just as the canonical transformation guarantees that a modified energyẼ h does exist, there might be a criterion for "compatibility" of vector fields with a canonical transformation. The next question would be whether such a "compatibility" holds only for the physical entropy or for all possible entropies (which would be the original degeneracy). Alternatively, as symplectic integrators can be obtained most easily from a variational principle [4,9,41,65], it might be worth looking at irreversible equations with a variational principle for GENERIC integrators.