{'fraction_non_alphanumeric': 0.05266093941760571, 'fraction_numerical': 0.017237532076313733, 'mean_word_length': 4.221963297407515, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Energy stable interior penalty discontinuous Galerkin finite element method for Cahn-Hilliard equation

An energy stable conservative method is developed for the Cahn--Hilliard (CH) equation with the degenerate mobility. The CH equation is discretized in space with the mass conserving symmetric interior penalty discontinuous Galerkin (SIPG) method. The resulting semi-discrete nonlinear system of ordinary differential equations are solved in time by the unconditionally energy stable average vector field (AVF) method. We prove that the AVF method preserves the energy decreasing property of the CH equation. Numerical results confirm the theoretical convergence rates and the performance of the proposed approach.


Introduction
In this work, we consider Cahn-Hilliard (CH) equation in mixed form as a system of second order differential equations in a bounded domain Ω ⊂ in Ω × (0, T ] (1a) in Ω × (0, T ] (1b) u(x, 0) = u 0 , in Ω (1c) with the homogenous Neumann boundary conditions ∂u ∂n = µ(u) ∂w ∂n = 0, on ∂Ω × (0, T ], or with periodic boundary conditions, where ∆ denotes the Laplace operator and µ(u) is the non-negative mobility function. The parameter provides a measure of the width of inter-facial layer, capturing the dominating effect of the reaction kinetics and represents the effective diffusivity. The CH equation is the most known model for phase separation and coarsening phenomena in a melted alloy [3]. It was also used as diffuse interface model for problems from fluid dynamics, material science and biology [4]. The variable u denotes the concentration of one of the species of the alloy and it is known as the phase state between materials.
In (2), the term δE(u)/δu stands for the variational derivative of the Ginzburg-Landau energy functional E given by with the potential function F (u) satisfying F (u) = f (u). In the literature two different types of potential functions are used: the convex double-well potential and the non-convex logarithmic potential [2,5] with 0 < θ ≤ θ c , where θ c is the transition temperature. For temperatures θ close to θ c , the logarithmic potential functional (5) is usually approximated by the quartic double-well potential (4). Other types of logarithmic potential functions can be found in [6,7,8]. Both potential functions (4) and (5) satisfy the monotonicity and the Lipschitz continuity conditions [9] ( for u 1 , u 2 ∈ R d with the constraints |u 1,2 | ≤ 1 and with C 1 , L f , L f ≥ 0. The two important properties of the CH equation are the conservation of mass under Neumann/periodic boundary conditions [1,2,4] Ω u(t)dx = Ω u 0 dx, and the decrease of total energy monotonically as t → ∞ In short, the CH equation describes a conservative gradient flow. The nonnegative mobility function µ(u) can be constant or degenerate [10], where the diffusion process is restricted to the interface zone, i.e. u ∈ [−1, 1] and it is zero outside. Commonly adapted versions of the degenerate mobility functions are µ(u) = βu(1 − u) [2,5] and µ(u) = β(1 − u 2 ) [6, 2, 4] for a positive β. From numerical point of view, it is desirable to develop numerical schemes which preserve the energy decreasing and mass conservative properties of the CH equation. Energy stability implies that the total energy of the fully discrete CH equation dissipates in time analogously to the continuous energy (3). The schemes that preserve the discrete versions of the continuous energy lead to approximate solutions which behave qualitatively similar to the continuous ones. Explicit methods are not suitable for time discretization of the CH equation because they are not energy stable and require very small time steps due to the stability restrictions.
In the literature, the CH equation with the constant mobility function has been discretized in space using finite differences [11], finite elements [7] and spectral methods [1]. Energy stable time discretization methods are based either on the convex splitting of the energy functional [2] or by adding stabilization terms to the energy functional F (u) [12]. On the other hand, the CH equation with degenerate mobility has been discretized by continuous finite elements [6,7], by local discontinuous Galerkin method [2,4,8], discontinuous Galerkin method with C 0 elements and with mixed finite elements [5], finite differences [13] and by spectral methods [14].
In this work, we use the mass conserving symmetric interior penalty discontinuous Galerkin (SIPG) finite element method [15,16] for the space discretization. Since the discontinuous Galerkin method uses piecewise polynomials, which are discontinuous at the interfaces, the discontinuous Galerkin (DG) approximation allows to capture the sharp gradients or singularities for small that affect the numerical solution locally. For the time integrator, it is well known that the first order backward Euler method is energy stable, i.e., the discrete energy decreases without any restriction of the step size ∆t for very stiff gradient systems with → 0 [17]. The only second order implicit energy stable method is the average vector field (AVF) method which preserves the energy decreasing property of the gradient systems and the systems with Lyapunov functions [18,17]. We prove that the energy decreasing property of the fully discrete CH equation with the degenerate mobility is preserved using SIPG in space and AVF method in time. Numerical convergence rates and numerical experiments for two dimensional problems with double well (4) and logarithmic (5) potential functions demonstrate the performance of the SIPG space discretization and the AVF time integrator.
The remainder of the paper is organized as follows. In the next section, Section 2, the SIPG discretization in space for the CH equation (1) is described and the semi-discrete system of ODEs in matrix-vector form is introduced. The fully discrete system using the AVF method as a time integrator is given in Section 3, and in Section 4, the energy decreasing property of the fully discrete system is proven. Numerical results are presented in Section 5 to demonstrate the accuracy of the numerical approach. The paper ends in Section 6 with some concluding remarks.

Semi-discrete formulations
In this section, we outline the semi-discretization of the CH equation (1) using symmetric interior penalty discontinuous Galerkin (SIPG), a type of discontinuous Galerkin (DG) methods, space discretization. To be being, the classical (continuous) weak formulation of the CH equation (1) reads as: for a.e. t ∈ (0, T ] where (·, ·) Ω denotes the usual L 2 -inner product over the domain Ω. It is well-known that under certain regularity assumptions, the system (6) has a unique solution in the space L 2 (0, T ; H 1 (Ω)).
In the sequel, we briefly describe the SIPG space discretization, and then, we introduce the matrix-vector form of the semi-discrete system of the CH equation (1) as a system of ODEs.

Discontinuous Galerkin discretization
In this work, we use as a space discretization method for the CH equation (1) the SIPG method [15,16] which is a type of discontinuous Galerkin finite elements method. Different from the classical finite elements methods, DG methods are suitable in the use of non-conforming grids, and requires lower regularity assumptions.
Let T h be a family of shape regular (triangular) elements The diameter of an element K and the length of an edge E are denoted by h K and h E , respectively. We set the space of discontinuous test and trial functions where P q (K) denotes the set of all polynomials on K ∈ T h of degree at most q. We split the set of all edges E h into the set E 0 h of interior edges and the set E ∂ h of boundary edges so that Let the edge E be a common edge for two elements K and K e . Since the functions u ∈ V h are discontinuous along the inter-element boundaries, there are two traces of u along E, denoted by u| E from inside K and u e | E from inside K e . Then, the jump and average of u across an interior edge E are defined, respectively, by where n K and n K e denote the outward unit normal vector to the boundary of the elements K and K e on the edge E, respectively. Similarly, for a vector function ∇u, the jump and average across an interior edge E are given by On a boundary edge E ⊂ ∂Ω, we set {u} = u and [u] = un, where n is the outward unit normal vector to the boundary. Then, in space SIPG discretized semi-discrete formulation of the CH equation (1) with a variable mobility function µ := µ(u) reads as: set where the bi-linear (in last two argument) form a h (κ; w, υ) is given by where the parameter σ is called the penalty parameter and it should be sufficiently large to ensure the stability of the SIPG discretization as described in [16] with a lower bound depending only on the polynomial degree [16]. In numerical experiments we take σ = 3q(q + 1).

Semi-discrete system in matrix-vector form
Let u h (t) and w h (t) be the time dependent solutions of the semi-discrete system (8). Then, u h (t) and w h (t) are of the form where ϕ m j are the basis functions spanning the space V h , and ξ m j and ζ m j are the unknown coefficients. The number n q stands for the local dimension depending on the basis polynomial order q and N is the number of (triangular) elements. By substituting the expansions (10) into the system (8) and choosing the test functions υ = ϕ k i , i = 1, . . . , n q , k = 1, . . . , N , we obtain the semi-linear systems of ordinary differential equations with the ordered unknown coefficient vectors and the basis functions In (11), M denotes the mass matrix with the entries M ij = (ϕ j , ϕ i ) Ω , 1 ≤ i, j ≤ n q × N , A µ and A are the stiffness matrices with the entries (

Fully discrete system
In this section, we give the fully discrete formulations of the CH equation (1) through the time integration of the semi-linear system of ODEs (11). To do this, we consider the uniform partition 0 = t 0 < t 1 < . . . < t J = T of the time interval [0, T ] with the uniform time step-size ∆t = t k − t k−1 , k = 1, 2, . . . , J.
In this work, as a time integrator we use the average vector field (AVF) method which is a structure preserving time integrator for gradient systems [18]. For a given system of ODEsẏ = −∇U (y), the AVF method is given by The AVF method possesses the energy decreasing property without restriction of the step sizes ∆t and is second order accurate in time. It represents a modification of the implicit mid-point rule and for quadratic potentials U (y), the AVF method reduces to the mid-point rule. Higher order variants of the AVF methods for Hamiltonian and Poisson systems with Gauss-Legendre collocation points are given in [17]. As Gauss-Legendre Runge-Kutta methods, the AVF method and higher order versions do not have damping property for very stiff systems, whereas for discontinuous Galerkin-Petrov methods and Radau II Runge-Kutta methods, the energy decreases monotonically without restriction of the step size ∆t and the Lipschitz constant for ∇U (y). However, they require the solution of coupled system of equations, which increases the computational cost for two and three dimensional CH equations, where efficient solution techniques are required [19]. Let for t = 0, u h,0 , w h,0 ∈ V h be the projections (orthogonal L 2 -projections) of the initial conditions u 0 , w 0 onto V h , and let η 0 = (ξ 0 , ζ 0 ) T be the corresponding initial coefficient vector each of components satisfying (10). At a specific time t = t n , we denote the coefficient vector of the approximate solutions (u h,n , w h,n ) T by η n = (ξ n , ζ n ) T . Then, the application of the AVF method (12) to the semi-discrete system (8) leads to solving for n = 0, 1, . . . , J − 1 the non-linear system of equations where we used the fact that the AVF method reduces the mid-point rule for linear terms. Further, the solution dependent mobility function µ(u h ) in the bilinear form A µ is computed explicitly as for continuous finite elements in [6], i.e., on each time interval (t n , t n+1 ] the mobility function is taken as µ = µ(u h,n ) where u h,n is the known approximate solution from the previous time step. The system (13) can be written as a residual equation R(η) = (R 1 (η), R 2 (η)) T = 0 with We solve the nonlinear system of equations (14) using the Newton's method. Starting with an initial guess η n+1 ) T , the k−th Newton iteration to solve the nonlinear system of equations (14) for the unknown vector η n+1 = (ξ n+1 , ζ n+1 ) T reads as until a user defined tolerance is satisfied. In (15), s = (s 1 , s 2 ) T is the increment, and J stands for the Jacobian matrix At each Newton iteration, the integral term in J b is approximated by the fourth order Gaussian quadrature formula.

Energy stability
It is expected that fully discrete energy stable schemes should preserve the discrete energy dissipation as their continuous parts, which leads to qualitatively better approximations. The continuous (in time) energy of the semidiscrete CH equation is given by [20] On the other hand, the discrete DG counterpart of the continuous energy (17) at a time t n = n∆t reads as In this section, we show that the AVF method applied to the semi-discrete system (8) is energy stable through the discrete energy (18). Applying the AVF method to the semi-discrete system (8), we obtain for any ν, ϑ ∈ V h (u n+1 − u n , ν) Ω + ∆t 2 a h (µ(u n ); w n+1 + w n , ν) = 0, Taking ν = (w n+1 + w n )/2 and ϑ = u n+1 − u n in (19), we obtain By using the identity (a + b, a − b) Ω = (a 2 − b 2 , 1) Ω and the bi-linearity of a h , we get Using the Taylor expansions of the terms F (u n ) and F (u n+1 ), and neglecting the higher order terms, we get We note that the bi-linear form a h satisfies a h (µ(u n ); w n+1 + w n , w n+1 + w n ) = µ(u n ) ∇(w n+1 + w n ) 2 Since all the terms in (23) are non-negative (see [16,Sec. 2.7.1] for positivity of edge integral term), we have a h (µ(u); w n+1 +w n , w n+1 +w n ) ≥ 0. Similarly, we have a h ( 2 ; u n+1 , u n+1 ) ≥ 0, a h ( 2 ; u n , u n ) ≥ 0. Using this identities, subtracting (20) from (21) and substituting (22), we obtain Hence, the AVF discretized scheme is energy stable through the discrete energy (18).

Numerical results
In this section, we present a set of numerical examples to confirm the accuracy, energy stability and mass conservation of our approach. All examples are considered in two dimensional spatial domain with constant and degenerate mobility function, double-well (4) and logarithmic (5) potential functions under Neumann and periodic boundary conditions.
We add a load vector to the system so that the exact solution above solves the system. The L 2 errors and the numerical order of accuracy at time t = 1 using the time step ∆t = 0.5∆x for linear and quadratic DG polynomials are given in Table 1. Both DG approximations satisfy the (k + 1)-th order of accuracy for P k elements.  [2] with the exact solution u(x, y, t) = e −2t sin(x) sin(y) in the domain Ω = [0, 2π] × [0, 2π] for t ∈ [0, 1] under periodic boundary conditions with the degenerate mobility function µ(u) = 1 − u 2 and doublewell potential (4). The effective diffusivity is taken as = 1. For first order DG polynomials we use time step ∆t = 0.0032π and for quadratic DG polynomials ∆t = 0.00032π. Order reduction is observed in Table 2 for linear and quadratic DG polynomials which is due to the nonlinearity of the degenerate mobility function.

Constant mobility and double-well potential under Neumann boundary conditions: spinodal decomposition and nucleation
We consider CH equation with constant mobility µ(u) = 1 and doublewell potential function (4) under homogenous Neumann boundary condition [2]. The computational domain is taken as Ω = [0, 1] × [0, 1] with = 10 −5 . This problem represents the two main separation mechanisms: spinodal decomposition and nucleation. Both mechanisms in the CH equation are defined by the initial condition u 0 (x) =ū + r whereū is a constant and r is a random number uniformly distributed on [0.005, 0.005]. Forū = 0, spinodal decomposition will be the governing separation mechanism (see Fig. 1). The mixture separates from a randomly perturbed homogeneous state (ū = 0) giving rise to a striped pattern of complicated topology that coarsens over time. If we let the simulation evolve, the stationary solution would be a fully separated flow with two rectangular patches. Whilē u = 0, the separation mechanism is nucleation (see Fig. 2). We show the solution forū = 0.4. In the nucleation mechanism, isolated nuclei come up from the mixture. Again, the spatial micro-structure of the mixture coarsens over time.
In both cases the discrete energy dissipates and the mass is conserved. Our results are similar to those in [21], where for spatial discretization local discontinuous Galerkin method and for time discretization implicit convex splitting have been used.

Degenerate mobility and logarithmic potential under Neumann boundary conditions
We consider 2D CH equation with degenerate mobility function µ(u) = u(1 − u) and the logarithmic potential function F (u) = 3000(u ln u + (1 − u) ln(1 − u)) + 9000u (1 − u) under homogenous Neumann boundary conditions [2] with the diffusion constant = 1 in the domain Ω = [−0.5, 0.5] × [−0.5, 0.5] for t ∈ [0, 0.2]. The initial condition is a random variation of uniform state u = 0.63 with a change no larger than 0.05. Fig. 3 shows the evolution of the concentration field. The two phases in the concentration evolution, the phase separation stage and the coarsening process stage can be seen clearly.

Conclusions
We have presented numerical results for the CH equation under periodic and homogenous Neumann boundary conditions using SIPG discretization in space and AVF method in time. The numerical energy is decreasing in all examples and the numerical results are in good agreement to the those in the literature.