Structure preserving reduced order modeling for gradient systems

Minimization of energy in gradient systems leads to formation of oscillatory and Turing patterns in reaction-diffusion systems. These patterns should be accurately computed using fine space and time meshes over long time horizons to reach the spatially inhomogeneous steady state. In this paper, a reduced order model (ROM) is developed which preserves the gradient dissipative structure. The coupled system of reaction-diffusion equations are discretized in space by the symmetric interior penalty discontinuous Galerkin (SIPG) method. The resulting system of ordinary differential equations (ODEs) are integrated in time by the average vector field (AVF) method, which preserves the energy dissipation of the gradient systems. The ROMs are constructed by the proper orthogonal decomposition (POD) with Galerkin projection. The nonlinear reaction terms are computed efficiently by discrete empirical interpolation method (DEIM). Preservation of the discrete energy of the FOMs and ROMs with POD-DEIM ensures the long term stability of the steady state solutions. Numerical simulations are performed for the gradient dissipative systems with two specific equations; real Ginzburg-Landau equation and Swift-Hohenberg equation. Numerical results demonstrate that the POD-DEIM reduced order solutions preserve well the energy dissipation over time and at the steady state.


Introduction
It is well known that energy minimization in gradient systems leads to formation of patterns. The energy driven pattern formation occurs in form of oscillatory and Turing patterns in gradient systems [1,2,3,4]. The longterm accurate computation of spatially heterogeneous steady state patterns require fine space and time meshes. Due to the implicit nature of time integrators used, a coupled fully nonlinear system has to be solved by the Newton-Raphson method at each time step with high accuracy. The discrete energy should be computed with high accuracy, when the nonlinear system is solved up to the machine precision. Therefore, the numerical solution of gradient systems for accurate resolution of patterns over long time integration might be very time consuming, especially for two and three dimensional (2D and 3D) problems.
The goal of this work is to develop an accurate and efficient reduced order model of 2D gradient systems by preserving the energy dissipation. During the last decade, model order reaction techniques have emerged as a powerful tool to reduce the cost of evaluating large systems of partial differential equations (PDEs). A low-dimensional linear reduced space is constructed that approximately represents the solution to the system of PDEs. Then, the solution can be obtained with a significantly reduced computational complexity by projection of the original system onto the reduced space (by Galerkin projection). Model reduction methods should produce low-dimensional surrogate models capable of mimicking the structural features of the underlying physics of the original system model. Reduced models that do not inherit key structural features and conservation properties of the original system may produce unphysical solutions.
The proper orthogonal decomposition (POD) is the most commonly used reduced order modeling technique in large-scale numerical simulations of complex systems. The stability of reduced order models over long-time integration and the structure preserving properties have been recently investigated in the context of Lagrangian systems [5,6], and for port-Hamiltonian systems [7,8]. For Hamiltonian and dissipative Hamiltonian systems, symplectic model reduction techniques are constructed in [9,10,11]. The average vector field (AVF) method is applied as energy preserving/dissipating time integrator to Hamiltonian and dissipative PDEs [12] like the nonlinear Schrödinger equation, Korteweg-de Vries equation, Sine-Gordon equation, Allen-Cahn equation, Ginzburg Landau equation and Maxwell equation. Applying AVF with reduced order modeling leads to stable and efficient solutions for energy dissipating systems like Allen-Cahn equation [13], FitzHugh-Nagumo equation [14] and energy conserving Hamiltonian systems [15,16].
There are a few works in the literature about the reduced order modeling of pattern formation in reaction-diffusion systems, like the complex Ginzburg-Landau equation [17,18] and FitzHugh Nagumo equation [14].
For the space discretization, we apply symmetric interior penalty discontinuous Galerkin (SIPG) method [19,20], and for time discretization, we utilize the energy preserving AVF method. The full order model based on SIPG-AVF discretization in space and time preserves the gradient structure of the Allen-Cahn equation [21], Cahn-Hilliard equation [22] and the skewgradient structure of the FitzHugh-Nagumo equation [14]. The POD is the most commonly used reduced order modeling technique in large-scale numerical simulations of complex systems. The POD-Galerkin approach involves an offline-online splitting strategy. In the offline stage, the high fidelity or truth solutions are generated by numerical simulations of the discretized high dimensional full order model (FOM). The POD is then applied to compute an optimal subspace to fit the high fidelity data. A reduced system is constructed by projecting the high-dimensional system to this subspace. In the online stage, the reduced system is solved in the low-dimensional subspace. The primary challenge in producing the low dimensional models of the high dimensional discretized PDEs is the efficient evaluation of the nonlinearities on the POD basis. Different POD Galerkin methods are developed to reduce the complexity of evaluating the nonlinear terms; gappy POD [23,24], missing point estimation [25,26], empirical and discrete empirical interpolation method (EIM/DEIM) [27,28]. Here, we apply the DEIM with the SIPG space discretization which is more efficient for evaluations of the nonlinear terms than the continuous FEM [14] in reduced order models. Moreover, efficient computation of the reduced order basis functions is an important aspect in reduced order modeling. The reduced order basis functions in POD-DEIM are computed in the offline stage using the singular value decomposition (SVD), which can be computationally demanding for large snapshot matrices. Here, we use the randomized singular value decomposition (rSVD) [29,30,31] as a fast and accurate alternative to the deterministic SVD to reduce the computational cost in the offline stage [32].
In the literature, there are some recent works concerning the ROMs based on POD which preserve energy dissipative structure of the model. For example, the study [16] offers mass and energy preserving FOM and ROM for two dimensional nonlinear Scrödinger equation (NLSE). To tackle the nonlinearity, DEIM is applied and it is shown that the discrete energy is preserved approximately for DEIM-ROM approximation. NLSE is solved by decomposing the real and complex parts of the solution which leads to a system of nonlinear PDEs. However, it is a Hamiltonian ODE. In another work [13], an energy stable ROM for a parameterized Allen-Cahn equation is proposed based on POD-greedy adaptive sampling of the snapshots in time and the greedy DEIM for the nonlinear term. It is reported that the energy decreases almost unconditionally for POD-DEIM approximation with large time stepsize. On the other hand, a skew-gradient system, namely FitzHugh-Nagumo equation, is investigated and a ROM-DEIM approximation is obtained and compared with FOM [14]. It is justified that the mini-maximizer property of the continuous energy holds for the full-discrete model. As different from the works mentioned above, in this current study, we investigate the structure preserving reduced order modeling for two typical gradient systems; Ginzburg-Landau equation and Swift-Hohenberg equation. Decomposing the solution of Ginzburg-Landau equation in real and complex parts, we obtain a system of parabolic-parabolic equations with cubic nonlinearities. On the other hand, the fourth order Swift-Hohenberg equation can be reduced to a system of parabolic-elliptic equations and the parabolic equation contains a cubic nonlinear term. After obtaining the systems, the full-discrete model is obtained by applying SIPG in space and AVF in time. Then, we derive the corresponding parabolic-parabolic and parabolic-elliptic ROMs with the use of DEIM for the nonlinearities. We prove the energy dissipation of the POD and POD-DEIM reduced order models and the contribution of the parabolicparabolic and parabolic-elliptic forms is reflected to the energy analysis. In detail, energy decreasing property of the gradient systems are proven under some condition on time step-size and this condition is related to both of the solutions of the parabolic-parabolic system, while the bound on time step-size depends only on the solution of the parabolic PDE for the parabolic-elliptic system. Moreover, we mention a relation between FOM and ROM so that the energy dissipation of FOM can be deduced easily. Furthermore, we compare DEIM and POD with respect to accuracy and efficiency; the POD-DEIM solutions have almost the same accuracy as the POD, but POD-DEIM is faster than POD. The steady state solutions of the gradient possesses Hamiltonian structure [2,3,4]. Our numerical results show that around the steady state, the energy of the FOMs and reduced order models (ROMs) remain constant. We drive upper bounds for the time meshes for the POD-DEIM reduced order models of the two gradient systems above. To the best of our knowledge, this is the first study justifying energy dissipation of POD-DEIM approximation for two typical example of gradient systems; Ginzburg-Landau equation and Swift-Hohenberg equation.
The remainder of the paper is organized as follows. The gradient systems are reviewed briefly in Section 2. Space-time discretization by SIPG and AVF method is given in Section 3. In Section 4, a small dimensional ROM by utilizing the POD method and also the POD-DEIM are presented. In Section 5, we investigate the gradient dissipative structures of reduced order RGL and SH systems based on DEIM approximation. The preservation of the dissipative structure, efficiency and accuracy of the ROMs are illustrated through the numerical simulations in Section 6. The paper ends with some conclusions.

Gradient systems
In the following section, we summarize the properties of gradient systems [2,3,33]. Let us consider an m-component reaction-diffusion system where w(x, t) = (w 1 (x, t), . . . , w m (x, t)) T and ∆ is the Laplace operator acting component-wise on w, t is the time variable, T > 0 is the final time and Ω ∈ R d (d = 1, 2, 3) is a bounded polygonal domain. We consider twodimensional reaction diffusion equations (d = 2) and set the spatial variable as x = (x, y). In (1), F ( w) is a real-valued potential functional, S is a nonnegative diagonal matrix, and Q is a non-degenerate symmetric matrix with Q 2 = I and satisfies the condition D T Q = QD, which ensures the nondegeneracy of QD. For the nonlinear vector f ( w) = (f 1 ( w), . . . , f m ( w)) T , the Jacobian matrix f w satisfies For the system (1), the corresponding energy functional is with ·, · R n denoting the standard Euclidean inner product on R n , and it holds The system (1) has a gradient structure when QS is a nonnegative symmetric matrix. In gradient systems, like the real Ginzburg-Landau equation and Swift-Hohenberg equation, (2) is the free energy which decreases monotonically, and local minimizers of the free energy correspond to the stable steady states. Below, we give two representative examples of the gradient systems, which are investigated in this paper.
Real Ginzburg-Landau (RGL) equation or Ginzburg-Pitaevski equation [34,35,36,2,3,37,38,39] is given by where the unknown solution ψ is complex-valued. Decomposing the solution ψ = u + iv in real and complex parts yields the coupled system of reactiondiffusion equations and it is of the form of the system (1) through the setting S = D = Q = I 2 , where I 2 is the two dimensional identity matrix, w = (u, v) T , and with the potential function Swift-Hohenberg (SH) equation [40,41,2,3,42,43,44,45] is given by Putting v = u + ∆u yields the coupled system In the sense of the system (1), the above system (4) is related to the setting The models (3) and (4) have cubic nonlinear reaction terms. Turing patterns are formed, when the critical wave number near the stationary solution is non-zero. In this case QD is indefinite. [2,3]. For RGL equation, QD is positive definite, so there exist no Turing patterns. However, there exist spatially periodic patterns around the stationary solutions [2,3]. Turing instability conditions and selection of Turing patterns for gradient systems are investigated in [3].
Steady state solutions of the gradient systems admit Hamiltonian structure [2,3,4]. The stationary system can be written in Hamiltonian form as with the skew-symmetric matrix J and the Hamiltonian given by The RGL equation (3) and SH equation (4) have dissipative structures for homogeneous Dirichlet and periodic boundary conditions [35,37,46,38]. Here, we consider periodic boundary conditions in both components for RGL and SH equations.

Discretization of the system (high-fidelity model)
In this section, we give a brief construction to the space-time discretization of the generic model (1), which we call high-fidelity or full order model (FOM). To do this, let us first consider the dimensionless scalar case of the generic model (1) as where Ω ⊂ R 2 is a bounded polygonal domain, and w 0 (x) is a prescribed initial condition. Later, we extend to the generic model (1) in a straightforward manner. In the sequel, we first give the semi-discretization of the system (5) by using SIPG method, which is a well-known member of discontinuous Galerkin (dG) methods, for space discretization, and then we present the full discretization using AVF time integrator.

Discontinuous Galerkin discretization in space
In the last twenty years, the dG methods gained an increasing importance for an efficient and accurate solution of PDEs. Although continuous finite elements methods (FEMs) require continuity of the solution along element interfaces, dG methods do not require continuity of the solution along the interfaces. In contrast to the stabilized continuous FEMs, dG methods produce stable discretization without the need for extra stabilization strategies and damp the unphysical oscillations. Due to the local structure of dG discretization, dG methods are parallelizable and adaptive meshing techniques can be implemented efficiently. The dG methods combine the best properties of the finite volume and continuous finite element methods.
To obtain the infinite dimensional variational formulation, as usual, we multiply the equation (5) by a test function φ(x) ∈ W , we take integral over the domain Ω, and we apply Green's formula to the second order term. Then, the continuous weak formulation of the problem (5) can be written as follows: for a.e. t ∈ (0, T ], find w := w(t, x) ∈ W such that for any φ : where ·, · := ·, · L 2 (Ω) denotes the L 2 -inner product on the space L 2 (Ω) of square integrable functionals over the domain Ω ⊂ R 2 , and together with the corresponding norm · := · L 2 (Ω) : Throughout this paper, unless otherwise stated, ·, · and · denote the standard L 2 inner product and norm, respectively, given above. The solution space W is the set as the space of functions in the Sobolev space H 1 (Ω) satisfying the periodicity on the boundary ∂Ω of the domain Ω as For the finite dimensional problem, we apply SIPG, a member of dG methods. Consider {T h } h as a family of shape-regular simplicial triangulations of Ω for h > 0. Each mesh T h consists of closed triangles such that Ω = K∈T h K holds. The diameter of an (triangular) element K and the length of an edge E are denoted by h K and h E , respectively. We split the set of all edges E h into the set E 0 h of interior edges and the set E p h of (periodic) boundary edge-pairs. An individual element of the set E p h is of the is the corresponding periodic edge-pair of E l with l > m. Let the edge E be a common edge for two elements K and K e . For a piecewise continuous scalar function w, there are two traces of w along E, denoted by w| E from inside K and w e | E from inside K e . Then, the jump and average of w across the edge E are defined by where n K (resp. n K e ) denotes the unit outward normal to ∂K (resp. ∂K e ). Similarly, for a piecewise continuous vector field ∇w, the jump and average across the edge E are given by where, for l > m, we use the fixed unit outward normal vector n l to the boundary edge E l ⊂ ∂Ω.
For continuous FEMs, the idea is to approximate the solutions w of (6) from a conforming, finite dimensional space W h ⊂ W . On the other hand, we point out that the space of solutions or test functions consist of piecewise discontinuous polynomials in dG methods. That is, no continuity constraints are explicitly imposed on the state and test functions across the element interfaces. As a consequence, weak formulations must include jump terms across interfaces, and typically penalty terms are added to control the jump terms. We define the space of solution and test functions by where P q (K) is the set of polynomials of degree at most q in K. Note that the space of discrete solutions and the space of test functions are identical due to the weak treatment of boundary conditions in dG methods. Note also that the space W h is a non-conforming space, i.e., Then, the SIPG weak formulation of the system (6) reads as: for a.e. t ∈ (0, T ], find where the bi-linear form a h (·, ·) and the non-linear form b h (·, ·) are given for where the bilinear form J p h (w h , φ) is related to the periodic boundary edges. Moreover, the parameter κ in the above formulation is called penalty parameter which should be sufficiently large to ensure the stability of the SIPG discretization with a lower bound depending only on the polynomial degree q, see, e.g., [20].
The SIPG discretized semi-discrete solution of (9) is given by where w i j (t) and ϕ i j (x) are the unknown coefficients and the basis functions in W h , respectively, for j = 1, 2, · · · , n q and i = 1, 2, · · · , n K . The number n K denotes the number of (triangular) elements in T h , and n q is the local dimension on each element with the identity n q = (q + 1)(q + 2)/2, where q is the degree of the polynomial order. Note that the degrees of freedom in dG methods are given by N := n K × n q , and throughout this paper we denote the dimension of the high-fidelity model or FOM by N . Inserting the expansions (10) into the system (9), we obtain the following dynamical system where w := w(t) ∈ R N is the unknown coefficient vector for the solution w h with the ordered entries w = (w 1 1 (t), · · · , w 1 nq (t), · · · , w 1 n K (t), · · · , w n K nq (t)), the matrix M ∈ R N ×N is the usual mass matrix, A ∈ R N ×N is the stiffness matrix corresponding to the bilinear form a h (·, ·), and b(w) ∈ R N is the vector corresponding to the non- where we have used the same ordering for basis functions as for the unknown coefficients

Time discretization
Energy stable time discretization methods preserve the dissipative structure of the numerical solution of gradient flow equations. The small values of the diffusion parameter lead to stiff systems after spatial discretization. Implicit/implicit-explicit methods are developed, since the explicit methods are not suitable for stiff systems and the fully implicit systems require solution of non-linear equations at each time step. In the semi-implicit schemes, the linear stiff part is treated implicitly and the nonlinear part is treated explicitly, so that at each time step a linear system of equations is solved. However, these methods do not preserve the energy dissipation of the system. Implicit Euler method and average vector field (AVF) method are energy stable time discretization techniques which are robust with small diffusion parameters. The AVF method is the only second order implicit energy stable method and it preserves energy decreasing property for the gradient systems.
We split the time interval [0, T ] into J equal length subintervals (t k−1 , t k ] as 0 = t 0 < t 1 < . . . < t J = T with the uniform step-size ∆t = t k − t k−1 , k = 1, 2, . . . , J. We set w n ≈ w(t n ) as the approximate solution vector at the time instance t = t n . For t = 0, let w h (0, x) ∈ W h be the projection (orthogonal L 2 -projection) of the initial condition w 0 (x) onto W h , and let w 0 be the corresponding initial coefficient vector satisfying the expansion (10). The AVF method for an arbitrary ODEẏ = g(y) is given by where y n ≈ y(t n ) for n = 1, . . . , J and y 0 = y(t 0 ). Then, applying the AVF scheme (12) to the semi-discrete system (11), the full discrete problem of (5) is read as: for n = 0, 1, . . . , J − 1, find w n+1 satisfying

Full discretization of the model
In this section, we extend the full discrete formulation of the dimensionless scalar model (5) to the generic model (1), and also we give the FOM formulations of the representative examples like RGL equation (3) and SH equation (4).
Let us consider the generic model (1) with m-component unknown solution w(x, t) = (w 1 (x, t), . . . , w m (x, t)) T together with its semi-discrete solu- (10), and we define the total coefficient vector as w(t) = (w 1,h (t) T , . . . , w m,h (t) T ) T ∈ R mN . Moreover, for each i = 1, . . . , m, let w i,n ≈ w i,h (t n ) denote the approximate coefficient vectors at the time instance t = t n , and approximate total coefficient vector is w n ≈ (w T 1,n , . . . , w T m,n ) T . Finally, we set w 0 = (w T 1,0 , . . . , w T m,0 ) T , where each w i,0 is the L 2 -projection of the initial conditions w i,0 (x)'s onto the dG space W h . Then, the FOM of the model (1) is given by where ⊗ denotes the Kronecker product, and b(w) = (b 1 (w), . . . , b m (w)) T ∈ R mN is the nonlinear vector whose entries are given by According to the general formulation (13), and the related settings of S, D, Q and potential functions F ( w) for the representative examples, their FOM formulations for the unknowns w := (u, v) T , their coefficients vectors w := (u, v) T , and the corresponding dG counter-part of the discrete energy functionals (the term D∇ w, Q∇ w R n in the energy formulation (2) corresponds to the dG bilinear form a h (·, ·) ) are expressed as FOM for SH equation (4) Note that in the SH equation (15), we do not require a nonlinear vector b 2 (u, v) since the function f 2 (u, v) for the corresponding system is linear in both u and v.

Reduced order model
Because the computation of the FOM (13) is time consuming, in this section, we construct a small dimensional ROM by utilizing the POD method [48]. In addition to POD, the nonlinear vectors in the ROM are computed efficiently by DEIM. The low-rank approximation is computed in three steps: computation of the numerical solutions of the original high-fidelity system; dimensionality-reduction of the snapshot matrices by SVD; Galerkin projection of the dynamics on the low-rank subspace. The first two steps are known as the offline stage, and the last one is the online stage. Offline stage is usually expensive and online step should be fast to run in real time.
As we did for dG formulation, we give, firstly, a brief construction procedure of the ROM for the scalar system (5) with the related dynamical system (11), and then we give the ROM formulations for our representative examples.

POD Galerkin projection
For the scalar system (5) with the N -dimensional FOM (11), the ROM of lower dimension k N is formed by the Galerkin projection of the system onto a k-dimensional reduced space resulting in the reduced solution w r h (t, x) ≈ w h (t, x), an approximation to the high-fidelity solution w h (t, x), and it satisfies the weak formulation The reduced solution w r h (t, x) is now of the form where ψ i (x)'s are the orthogonal (in L 2 -sense) reduced basis functions spanning the reduced space W r h , and w r i (t)'s are the coefficients of the reduced solution, which we call reduced coefficients. From the reduced coefficients w r i (t), we set the solution of the reduced system as w r (t) := (w r 1 (t), . . . , w r k (t)) T . Note that the reduced space W r h ⊂ W h is a subset of the dG space W h ; hence the reduced basis functions are linear combination of the dG basis functions {ϕ j } N j=1 : where Ψ j,i 's are the coefficients of the i-th reduced basis function ψ i (x). Then, using the column vectors Ψ i = (Ψ 1,i , . . . , Ψ N,i ) T of the coefficients, we construct the following matrix of POD modes To obtain the reduced basis functions {ψ i (x)} k i=1 , we need to solve the minimization problem [48] min where δ ij is the Kronecker delta. The above minimization problem is equivalent to the eigenvalue problem [48] where W := [w 1 , . . . , w J ] ∈ R N ×J is the snapshot matrix whose n-th column vector w n is the solution vector of the FOM at time t n . Then, the matrix Ψ of POD modes can be computed through the SVD of the snapshot matrix W.
In addition, between the coefficient vectors w(t) of the high-fidelity solutions and the coefficient vectors w r (t) of reduced solutions, we have the relation from where we can find the initial reduced vector as w r (0) = Ψ T M w 0 . For a more detailed description, we refer to the study [14].
Projecting the system (16) onto the reduced space W r h , and by the use of above formulations, we finally obtain the following ROM as a dynamical system: with the reduced stiffness matrix and the reduced nonlinear vector Like the FOM (11), the ROM (20) is solved in time by the AVF method.

Approximation of the nonlinearities
Although the reduced model (20) is of small dimension, the computation of the nonlinear vector b r (w r ) still depends on the dimension N of the FOM. In this section, we explain DEIM [28] to reduce the computational complexity due to the nonlinear vector in the ROM (20).
The DEIM aims to find an approximation to the nonlinear vector b(Ψw r (t)) which is the full dimensional nonlinear part of the reduced nonlinear vector b r (w r ) = Ψ T b(Ψw r ) for a.e. t ∈ (0, T ]. It is followed by the projection onto a subspace of the column space of the (nonlinear) snapshot matrix B := [b(w 1 ), . . . , b(w J )] ∈ R N ×J whose columns are the nonlinear vectors at discrete times t n for n = 1, . . . , J. Then, a DEIM basis {Q i } m i=1 of dimension m N is constructed by applying the POD to the matrix B, and the following approximation is used: where Q := [Q 1 . . . Q m ] ∈ R N ×m is the DEIM basis matrix, and c(t) is the unknown coefficients vector to be determined. Since the system (21) is overdetermined, a projection matrix P = [e p 1 , . . . , e pm ] ∈ R N ×m with e p i is the i-th column of the identity matrix I ∈ R N ×N is computed. Then, the approximate nonlinearity is derived aŝ b(Ψw r (t)) = Q(P T Q) −1 P T b(Ψw r (t)), (22) and then the reduced model (20) can be rewritten as: where the matrix B := Ψ T Q(P T Q) −1 is computed once in the off-line stage, and the reduced nonlinear vector b r deim (w r ) := P T b(Ψw r ) requires only m N integral evaluations. On the other hand, the computation of the Jacobian of the nonlinear vector requires N × n p integral evaluations without DEIM, but it is only m × n p with DEIM approximation.
For the details of the computation of the reduced non-linear vectors, we refer to the greedy DEIM algorithm [28]. For continuous finite element and finite volume discretizations, the number of flops for the computation of bilinear form and nonlinear term depends on the maximum number of neighbor cells [49]. In the case of dG discretization, due to its local nature, it depends only on the number of nodes in the local cells. For instance, in the case of SIPG with linear elements (n q = 3), it contributes only 3 nonzero integrals for each degree of freedom, hence 3 integrals have to be computed on a single triangular element [50], whereas in the case of continuous finite elements, integral computations on 6 neighbor cells are needed [51]. Since the AVF method is an implicit time integrator, at each time step, a non-linear system of equations has to be solved by Newton's method. The reduced Jacobian has a block diagonal structure for the SIPG discretization, which is easily invertible [50], and requires only O(n q N ) operations with DEIM.

ROM systems for the representative examples
In this section, we give the ROM formulations related to the FOM systems (14) and (15) of RGL and SH equations, respectively. In all these two systems, the unknown solution is w = (u, v) T having two components. So, for a.e. t ∈ (0, T ], let us consider the dG (high-fidelity) solutions u h (t, x) ∈ W h and v h (t, x) ∈ W h with the related FOM solutions (coefficient vectors of the solutions) u(t) ∈ R N and v(t) ∈ R N . Then, denote the corresponding reduced coefficient vectors by u r (t) ∈ R ku and v r (t) ∈ R kv , which satisfy where Ψ u ∈ R N ×ku and Ψ v ∈ R N ×kv are the matrices of POD modes related to the snapshot matrices of coefficient vectors of the high-fidelity solutions u h (t, x) and v h (t, x), respectively. Thus, the columns of the matrices of POD modes Ψ u and Ψ v are the coefficient vectors of the reduced basis functions ψ u (x) and ψ v (x) for the high-fidelity solutions u h (t, x) and v h (t, x), respectively, and they are the solutions of the eigenvalue problems (19) with W ∈ {U, V}, Ψ i ∈ {Ψ u,i , Ψ v,i }, and k ∈ {k u , k v }. Moreover, let Q 1 ∈ R N ×m 1 and Q 2 ∈ R N ×m 2 denote the DEIM basis of dimension m 1 and m 2 , respectively, of the following (nonlinear) snapshot matrices: with the corresponding DEIM projection matrices P 1 and P 2 . Under these assumptions, the ROM with DEIM formulations for the RGL and SH equations are given as the following systems:

Energy analysis of ROM with DEIM approximation
In this section, we investigate the gradient dissipative structures of reduced order models for RGL and SH systems based on DEIM approximation. We examine the difference of energy functionals at two successive time steps, , and we find some upper bounds for it using the corresponding full discrete ROMs. The key point in the analysis is to justify the equivalence of the integrals of f 1 (u, v) and f 2 (u, v), which appear in the full discrete weak form due to AVF method, and the potential F (u, v). Moreover, we include the contribution of DEIM by adding and subtracting the nonlinear terms to the difference of energy functionals written in terms of the weak forms, since DEIM offers an approximation for the nonlinear terms instead of exactly evaluating them. Then, we can examine the effect of DEIM on dissipative structures. On the other hand, the contribution of the parabolic-parabolic and parabolic-elliptic forms is reflected to the energy analysis. In detail, energy decreasing property of the gradient systems are proven under some condition on time step-size and this condition is related to both of the solutions of the parabolic-parabolic system, while the bound on time step-size depends only on the solution of the parabolic PDE for the parabolic-elliptic system. Moreover, we note that in the preceding analysis, the AVF method applied to linear terms coincides with the mid-point rule.

Real Ginzburg-Landau system
We write the corresponding reduced full discrete (AVF method is applied) variational formulation associated with the FOM (14) of the RGL equation as where, for analysis and computational purposes, we exclude the linear terms µu and µv from the nonlinearities, and we have set Choosing the test functions as φ = u r n+1 − u r n and η = v r n+1 − v r n from the reduced space W r h , and using the identity (a + b)(a − b) = a 2 − b 2 on dG bilinear forms and linear terms on the right hand sides, we get Now, we use the change of variables to obtain Through the use of (28), the side-by-side summation of the equations in (27) gives rise to where we used the fact that Hence, the identity (29) implies that E h [u r n , v r n ]−E h [u r n+1 , v r n+1 ] > 0, and that without the use of DEIM approximation for the nonlinearities, the discrete energy of the ROM for RGL equation decreases as the time progresses.
Remark 5.1. The identity (29) implies that not only the ROM without DEIM approximation preserves the discrete energy decreasing property exactly, but also FOM (14) of RGL equation preserves exactly the decreasing structure of the discrete energy. To see this, it is enough just to replace the reduced solution space W r h to the full solution space W h , and hence u r h and v r h are replaced with u h and v h , respectively, in the above analysis.
In the case of the ROM (24) of RGL equation with the use of DEIM approximation for the nonlinearities, we should proceed further. Before investigating the energy functional, we particularly consider the integral term in (27) containing the nonlinearityf 1 (ũ r h ,ṽ r h ) to obtain the following relation 1 0 Ωf The relation for the integral term containing the nonlinearityf 2 (ũ r h ,ṽ r h ) can be found in a similar manner. In order to catch the contribution of DEIM approximations, we need to replace the nonlinear vectorsb 1 andb 2 with their DEIM approximationsb 1 andb 2 , respectively, given in (21). Now, we apply the replacements of DEIM approximations, we add and subtract the terms Ωf 2 (ũ r h ,ṽ r h )(v r n+1 − v r n ) dx dξ to the right-hand side of (27), and we obtain the following relation [13,Sec.3 where the last two rows on the right hand side are additional because of the DEIM contribution. In particular, with the use of the integral mean theorem and the a priori estimates in [28,52,51] for the DEIM approximation, we further obtain for someû r n between u r n and u r n+1 , where the matrix R is the Cholesky factor of the mass matrix M , and · 2 denotes the usual Euclidean inner product on R N . Applying the same idea forb 2 (Ψ uû r (t), Ψ vv (t)), combining (31) and (32), we finally obtain We deduce that ROM with DEIM formulation preserves the energy decreasing property, i.e., (34) By construction, the matrices Q 1 and Q 2 are orthonormal, and the terms (P T 1 Q 1 ) −1 2 and (P T 2 Q 2 ) −1 2 are of moderate magnitude. They differ in the numerical tests between 10 − 30, and 30 − 60, respectively. Hence, the upper bound for the time step-size ∆t in the right hand side of (34) is a sufficiently large number so that the discrete energy for the ROM (24) of RGL equation with the use of DEIM approximation decreases almost unconditionally for large time step-sizes.

Swift-Hohenberg system
The AVF applied reduced variational formulation related to the FOM (15) of the SH equation is given by where, again for computational purposes as we did for RGL equation, we have set the nonlinearity f (u) = µu − u 3 , andũ r = ξu r n+1 + (1 − ξ)u r n . Choosing φ = u r n+1 − u r n and η = −(v r n+1 − v r n ) in (35), summing up the equations side-by-side, using the bi-linearity of a h (·, ·), and from the identity With the change of variables z :=ũ r = ξu r n+1 + (1 − ξ)u r n , dz = (u r n+1 − u r n )dξ, we obtain that Combining (36) and (37), we obtain that Clearly, the ROM without DEIM approximation, and also the FOM for SH equation (Remark 5.1), preserves the energy decreasing property exactly. For the contribution of the DEIM approximation, we follow the similar steps as we did for the RGL equation: we replace the nonlinear vector with its DEIM approximation, we add and subtract the term 1 0 f (ũ r h ), u r n+1 − u r n dξ to the right-hand side of (37) to reach at for someû n r between u n r and u n+1 r , and with the DEIM approximation b(Ψw r (t)) ≈b(Ψw r (t)). From (38), we find out that the discrete energy decreasing property, i.e., E h [u r n+1 , v r n+1 ] ≤ E h [u r n , v r n ], holds for ROM (25) of SH equation with the DEIM approximation under the following condition on the step-size which allows the energy decreasing property almost unconditionally, by the same reasoning we have stated for the RGL equation.

Numerical Results
In this section, we present some numerical results for the two gradient systems. In all simulations, we use linear dG basis functions on a uniform 32 × 32 rectangular grid with 2048 triangular elements in the rectangular domain Ω = [a, b] 2 . The number of POD and DEIM modes are determined by a user-defined threshold such that "relative information content" criterion satisfies that RIC(k) ≥ 1 − ε. In (40) σ i is the i-th singular value, and d z is the rank of the snapshot matrix Z. When ε = 0, the snapshot data is recovered perfectly up to numerical errors. In the numerical simulations, we took for the DEIM modes a smaller value of threshold than for the POD so that the POD and DEIM errors have the same level of accuracy. For the computation of the SVD, we use the randomized SVD (rSVD) which is computationally much faster then the standard SVD. For the details of rSVD algorithm we refer to the studies [53,16]. In the numerical experiments, we have used the same time-step size ∆t = 0.01 for the FOMs and ROMs. Only one Newton iteration is needed in solving nonlinear equations of the FOMs and ROMs. All simulations are performed on a Windows 10 machine with Intel Core i7, 2.5 GHz and 8 GB using MATLAB R2014. CPU time is measured in seconds. The numerical simulations are stopped until the system reaches a spatially inhomogeneous steady-state satisfying, for a user defined tolerance tol, the following criteria The errors between the FOM and ROM solutions are given in L 2 -L 2 norm (Frobenious norm) over the discrete space and time, whereas the energy errors are computed in L 2 norm over the time interval. We present the FOM and ROM solutions for the GL equation only for the u component, because the patterns for v are similar to u.
The singular values are plotted in semi-logarithmic form in Fig. 1, with the index of singular values on the horizontal axis. The singular values of the snapshot matrices and of the nonlinear terms decay linearly in Fig 1  . The threshold for the selection of POD and DEIM modes are chosen as = 10 −4 and = 10 −6 , respectively. In Fig. 2 the FOM solutions for 4 POD and 5 DEIM modes are shown at the steady state. For a small number of POD and DEIM modes, the FOM and ROM are very close in Fig. 2 and the L 2 -L 2 errors, given in Table 1, are relatively small. The FOM energy decays rapidly for a short time, and around T = 22 it remains constant until the steady-state in Fig. 3. The energy is very well preserved by the ROMs over the whole time integration as shown in Fig. 3 and Table 1. Using the DEIM modes, the computational cost is reduced about a factor of 7 over the POD solutions, see Table 2.
Decay of the singular values shown in Fig. 4 is similar to the RGL equation in Fig. 1, and the energy dissipates like RGL equation as presented in Fig. 6. However, the energy decay is not well preserved for the POD and DEIM. Because the DEIM energy errors are larger then POD energy errors, the reduced order patterns obtained by DEIM are less accurate than those obtained by POD, see Table 3 and Fig. 5. Because more DEIM modes are required than the POD modes, the speed-up of the DEIM is lower in Table  4 than for the RGL equation in Table 2.

Conclusions
In this paper, we have developed a structure preserving reduced order modeling for the two typical gradient systems. Numerical tests on twodimensional real Ginzburg-Landau and Swift-Hohenberg equations have illustrated the computational efficiency and the accuracy of the structure preserving reduced order models. The accuracy of the patterns and efficiency of the reduced order model depend on how well the energy evolution of the full order model is preserved by POD and DEIM. Because the singular values of the snapshot matrices and nonlinear terms decay slowly, relatively more POD and DEIM modes required. The slow decay of the singular values is a characteristic for advection dominated and transport problems with wave type solutions. In the recent years, new methodologies are developed to tackle this problem. Nevertheless, the time-dependent problems dominated by transport or propagation phenomena like the pattern formation remain a challenge for reduced order modeling. In a future study, we aim to apply these methodologies to reaction-diffusion systems with pattern formation.