a mixed discontinuous galerkin, convex splitting

DISCRETE AND CONTINUOUS DYNAMICAL SYSTEMS SERIES B Volume 18, Number 9, November 2013

doi:10.3934/dcdsb.2013.18.2211 pp. 2211–2238

A MIXED DISCONTINUOUS GALERKIN, CONVEX SPLITTING SCHEME FOR A MODIFIED CAHN-HILLIARD EQUATION AND AN EFFICIENT NONLINEAR MULTIGRID SOLVER

Andreas C. Aristotelous Statistical and Applied Mathematical Sciences Institute (SAMSI) Research Triangle Park, NC 27709, USA

Ohannes Karakashian and Steven M. Wise Department of Mathematics University of Tennessee Knoxville, TN 37996, USA

(Communicated by Jerry Bona) Abstract. In this paper we devise and analyze a mixed discontinuous Galerkin finite element method for a modified Cahn-Hilliard equation that models phase separation in diblock copolymer melts. The time discretization is based on a convex splitting of the energy of the equation. We prove that our scheme is unconditionally energy stable with respect to a spatially discrete analogue of the continuous free energy of the system, unconditionally uniquely solvable, and convergent in the natural energy norm with optimal rates. We describe an efficient nonlinear multigrid solver for advancing our semi-implicit scheme in time and conclude the paper with numerical tests confirming the predicted convergence rates and suggesting the near-optimal complexity of the solver.

1. Introduction. Let Ω ⊂ Rd , d = 2, 3, be an open polygonal or polyhedral domain. We consider the following modified Cahn-Hilliard problem with natural boundary conditions [10, 9, 25]: ∂t u = ∆w − θ (u − u0 ) , 1 3 u − u − ε∆u, w= ε ∂n u = ∂n w = 0,

in ΩT := Ω × (0, T ),

(1a)

in ΩT ,

(1b)

on ∂Ω × (0, T ),

(1c)

on Ω × {t = 0} , (1d) R −1 where ε > 0 and θ ∈ R are parameters, and u0 = |Ω| u0 (x) dx. The parameter ε is called the interfacial thickness (or just interface) parameter. When θ = 0, the classical Cahn-Hilliard equation is recovered [8, 14, 15]. Equations (1a) – (1b) have been proposed as a model for phase separation of diblock copolymer melts [9, 10, 25] and possess solutions that can be significantly different from those of the classical Cahn-Hilliard equation. See, for example, Figs. 1 and 5 and the works cited. Herein u = u0 ,

2010 Mathematics Subject Classification. Primary: 65M60, 65M12; Secondary: 65M55. Key words and phrases. Cahn-Hilliard equation, diblock copolymer, discontinuous Galerkin, symmetric interior penalty, convex splitting, energy stability, convergence, nonlinear multigrid. The second author is supported by NSF-DMS 0811314. The third author is supported by NSF-DMS 1115390.

2211

2212

A. C. ARISTOTELOUS, O. KARAKASHIAN AND S. M. WISE

t = 0.002, θ = 15000.0

t = 0.002, θ = 5000.0

t = 0.002, θ = 0.0

t = 0.01, θ = 15000.0

t = 0.01, θ = 5000.0

t = 0.01, θ = 0.0

t = 0.02, θ = 15000.0

t = 0.02, θ = 5000.0

t = 0.02, θ = 0.0

Figure 1. Approximate solutions to the modified Cahn-Hilliard equation (1a) – (1b) with three different values of θ. In column 1, θ = 15000.0; in column 2, θ = 5000.0; and in column 3, θ = 0.0. For all three computations u0 = −0.1, ε = 0.02, and Ω = (0, 8)2 . Observe that for large values of θ the coarsening process is significantly inhibited. The maximum (umax ) and minimum (umin ) values of u are also altered upon changing θ. In column 1, umax ≈ 0.75, umin ≈ −0.75; in column 2, umax ≈ 0.95, umin ≈ −0.95; and in column 3, umax ≈ 1, umin ≈ −1.

we call u the phase variable: u = 1 represents pure polymer A, and u = −1, pure polymer B. Physically realizable copolymer mixtures satisfy −1 ≤ u ≤ 1. The variable w is called the chemical potential, though it is not actually a chemical potential in the strictest sense, since it is not the variation of the energy that is given below. A weak formulation of (1a) – (1d) may be written as follows: find (u, w) such 1 4 ∞ 2 −1 that u ∈ L∞ 0, T ; H (Ω) ∩ L (0, T ; L (Ω)), ∂t u ∈ L 0, T ; H (Ω) , and w ∈ 2 1 L 0, T ; H (Ω) , and there hold for almost all t ∈ (0, T ) h∂t u, χi + ∇w, ∇χ + θ (u − u0 , χ) = 0

∀χ ∈ H 1 (Ω),

(2a)

A MIXED DG CS SCHEME FOR A MODIFIED CH EQUATION

2213

1 3 u − u, ψ = 0 ∀ψ ∈ H 1 (Ω), (2b) ε with the “compatible” initial data 2 u( · , 0) = u0 ∈ HN (Ω) := v ∈ H 2 (Ω) ∂n v = 0 on ∂Ω . (3) ∗ Here we use the notation H −1 (Ω) = H 1 (Ω) , and h · , · i is the duality paring between H −1 (Ω) and H 1 (Ω). The system (2a) – (2b) is mass conservative: for almost every t ∈ [0, T ], (u( · , t) − u0 , 1) = 0. Observe that the homogeneous Neumann boundary conditions are natural in this mixed weak formulation of the problem. The existence of weak solutions is a straightforward exercise using the compactness/energy method. See, for example, [18]. To define the energy for this system we need a norm on a subspace of H −1 (Ω). With L20 (Ω) denoting those function in L2 (Ω) with zero mean, we set ˚1 (Ω) = H 1 (Ω) ∩ L20 (Ω), H ˚−1 (Ω) := v ∈ H −1 (Ω) hv, 1i = 0 , (4) H (w, ψ) − ε (∇u, ∇ψ) −

˚−1 (Ω) → H ˚1 (Ω) via the following variational We define a linear operator T : H −1 1 ˚ ˚ problem: given ζ ∈ H (Ω), find T(ζ) ∈ H (Ω) such that ˚1 (Ω). (∇T(ζ), ∇χ) = hζ, χi ∀χ ∈ H (5) T is well-defined, as is guaranteed by the Riesz Representation Theorem. The following facts are easily established. ˚−1 (Ω) and, for such functions, set Lemma 1.1. Let ζ, ξ ∈ H (ζ, ξ)H −1 := (∇T(ζ), ∇T(ξ)) = (ζ, T(ξ)) = (T(ζ), ξ) .

(6)

˚−1

( · , · )H −1 defines an inner product on H (Ω), and the induced norm is equivalent to the operator norm: q hζ, χi . (7) kζkH −1 := (ζ, ζ)H −1 = sup k∇χk ˚1 L2 06=χ∈H ˚−1 (Ω), |hζ, χi| ≤ kζk −1 k∇χk 2 . Consequently, for all χ ∈ H 1 (Ω) and all ζ ∈ H H L 2 Furthermore, for all ζ ∈ L0 (Ω), we have the Poincare type inequality: for some C > 0, kζkH −1 ≤ C kζkL2 . Consider the energy 1 1 |Ω| ε θ 4 2 2 2 E(u) := kukL4 − kukL2 + + k∇ukL2 + ku − u0 kH −1 , (8) 4ε 2ε 4ε 2 2 which is defined for all u ∈ A := u ∈ H 1 (Ω) (u − u0 , 1) = 0 . (The affine set A is called the admissible set.) Clearly, if θ ≥ 0, then E(u) ≥ 0 for all u ∈ A. For any θ ∈ R, ε > 0, and u ∈ A, there exist positive constants M1 = M1 (ε, θ) and M2 = M2 (ε, θ) such that 2

0 ≤ M1 kukH 1 ≤ E(u) + M2 ,

(9)

i.e., the energy is coercive. Thus, one can prove that minimizers exist by the standard methods. It is straightforward to show that weak solutions of (2a) – (2b) dissipate the energy (8). In other words, (1a) – (1d) is a conserved gradient flow with respect to the energy (8). Precisely, for any t ∈ [0, T ], we have the energy law Z t 2 E(u( · , t)) + k∂t u( · , s)kH −1 ds = E(u0 ). (10) 0

2214

A. C. ARISTOTELOUS, O. KARAKASHIAN AND S. M. WISE 2

When θ ≥ 0, the energy term θ2 ku − u0 kH −1 in (8) is convex and stabilizing. The simulation results in Fig. 1 reflect this case. Physically, when θ ≥ 0 the term 2 θ 2 ku − u0 kH −1 tends to suppress the coarsening process. Otherwise, if θ ≤ 0 the term is concave and destabilizing. In this case (not shown), the process of phase separation should be enhanced. In the general case, we can write θ = θc − θe , where θc , θe ≥ 0. The energy then has the (non-unique) convex splitting E = Ec − Ee , where |Ω| ε θc 1 4 2 2 kukL4 + + k∇ukL2 + ku − u0 kH −1 , 4ε 4ε 2 2 (11) 1 θe 2 2 Ee (u) = kukL2 + ku − u0 kH −1 . 2ε 2 Ec and Ee are clearly convex. The convexity splitting scheme that we introduce later will respect this convex decomposition of the energy in order to guarantee stability and solvability properties. To make the presentation more manageable, we only consider the case that θ ≥ 0, i.e., θe = 0. However, the more general case may be of some theoretical and physical interest and is easily handled in the framework we develop. The goal of this paper is to construct a fully discrete, unconditionally energy stable, unconditionally uniquely solvable, and convergent discontinuous Galerkin finite element (DGFE) scheme for the modified Cahn-Hilliard equation (1a) – (1b). Energy stability means that the discrete solutions dissipate some suitable energy in time, in analogy with the PDE model. In other words, the scheme will preserve a discrete version of the continuous energy dissipation law (10). It is generally accepted that schemes that preserve discrete analogues of the energy dissipation law and/or other invariants of the PDE lead to approximations that behave qualitatively more like the true PDE solutions. In this work we utilize a mixed formulation of our problem and use a symmetric interior penalty discontinuous Galerkin finite element (SIP-DGFE) method for the spatial discretization [3]. For the time discretization we use the convex-splitting (CS) approach popularized by David Eyre [16] and used in several other settings [12, 18, 29, 33]. In the CS framework, one treats the contribution from the convex part implicitly and the contribution from the concave part explicitly in the scheme. This treatment automatically confers to the schemes energy stability and unique solvability, both properties being unconditional with respect to the time and space step sizes. Pioneering work in primitive variable discontinuous Galerkin methods for fourthorder elliptic problems was done by Babuˇska and Zl´amal [4] and Baker [5]. One of the principal motivations for these efforts was the possibility to avoid using C 1 spatial elements to approximate solutions to higher-order equations. DG methods for second-order elliptic and parabolic equations were developed subsequently in [2, 13, 26, 31]. While there has been extensive work for Cahn-Hilliard-type equations using the standard continuous galerkin finite element setting, only a few papers have investigated the use of discontinuous galerkin methods in this setting. The Cahn-Hilliard scheme of Feng and Karakashian [17] is of primitive variable, SIPDGFE type, rather than a mixed scheme like we propose and examine here. While the paper by Feng and Karakashian [17] does allow for a dynamic mesh, their scheme may not be energy stable or unconditionally solvable. They are able to prove convergence in the dynamic mesh setting, which remains an important result. The scheme in [30] is a primitive variable, C 0 -DGFE scheme, whereas the one here is totally discontinuous and of mixed type. The DGFE scheme for the CH equation Ec (u) =


2215

proposed in [34] is of LDG type, meaning (essentially) that there are variables for u and each of its three spatial derivatives. It does appear that the scheme in [34] is at least conditionally energy stable. The authors in [24] considered and analyzed a mixed fully implicit scheme for the Cahn-Hilliard equation with convection (i.e, with an added prescribed, convection velocity) using a backward Euler discretization in time. They proved convergence of their scheme, but only conditional solvability and conditional a priori stability. To our knowledge there has been no numerical analysis, in either the standard continuous Galerkin or discontinuous Galerkin finite element settings, for the the modified Cahn-Hilliard equation that we consider here. The SIP-DGFE framework that we employ for the modified Cahn-Hilliard equation (1a) – (1b) is of mixed type, involving the two natural variables of the problem u and w. Here we include the θ “growth” term that was not considered previously, an addition that makes the analysis and the results significantly different in the current presentation. Specifically, analysis of this term requires the introduction of a broken negative norm. Furthermore, unlike the results in [24], the energy stability and solvability statements we prove are completely unconditional with respect to the time and space step sizes since we use a CS framework. In fact, all of our a priori stability estimates hold completely independently of the time and space step sizes that are used. We also include a detailed proof of the unconditional unique solvability of our method, showing that the scheme is equivalent to a convex minimization problem, something that was not done in the previous works mentioned. In addition, we address the (non-trivial) solution of the algebraic system of equations that result from our approximation scheme using a novel, efficient nonlinear multigrid solver. Specifically, we show with a reasonable battery of tests that the iterative multigrid solver enjoys nearly optimal convergence to the solution of the nonlinear algebraic system of equations. We view these features collectively as a meaningful advance of the state-of-the-art for DGFE methods for the Cahn-Hilliard family of equations. Our results here do not hold if the meshes are locally and dynamically refined or coarsened in time in an arbitrary way. Also, we do not prove convergence of the multigrid solver. These deficiencies represent important lines of work that we plan to take up in the near future. For now, the interested reader can see [1]. The paper is organized as follows. In Section 2 we briefly review some notation and tools from the SIP-DGFE setting. In Section 3 we define the scheme and establish the unconditional solvability and stability of the scheme. In Section 4 we prove that the scheme is convergent under suitable regularity assumptions for the PDE solution. In Section 5 we describe a nonlinear multigrid solver for our scheme. In Section 6 we present the results of numerical tests that corroborate the convergence and energy stabilities of the scheme. Therein we also give the results of some tests that demonstrate the potential efficiency for our nonlinear multigrid solver in that section. 2. Discontinuous Galerkin preliminaries. Let Th = {K} be a family of starlike partitions of the domain Ω ⊂ Rd (d = 2 or 3), parametrized by 0 < h < 1, where h = maxK∈Th hK and hK = diam(K). We will assume that Th satisfies the following conditions: T1. The elements (cells) of Th satisfy a minimum angle condition, i.e, there exists a constant 0 < C0 ≤ 1, independent of h, such that C0 ρK ≤ hK , where ρK is the diameter of the inscribed ball to K.

2216


T2. Th is locally quasi-uniform, i.e., there exist constants 0 < C1 ≤ C2 , independent of h such that, for any two adjacent cells K, K 0 ∈ Th , C1 hK ≤ hK 0 ≤ C2 hK . (Two cells K, K 0 ∈ Th are called adjacent if and only if measd−1 (∂K ∩ ∂K 0 ) > 0.) T3. Th is globally quasi-uniform, i.e., there exist constants 0 < C3 ≤ C4 , independent of h such that, for any two (not necessarily adjacent) cells K, K 0 ∈ Th , C3 hK ≤ hK 0 ≤ C4 hK . We always enforce conditions T1 and T2. The third statement is stronger and more restrictive than the second, and we will point out when condition T3 is specifically needed. The DGFE framework involves functions that are discontinuous across interelement boundaries. This motivates the use of so-called “broken” Sobolev spaces: Y H m (Th ) := H m (K) = v ∈ L2 (Ω) v|K ∈ H m (K), ∀K ∈ Th . K∈Th

It is important to note that the members of H m (Th ) are not functions in the proper sense, since they can be multivalued on the inter-element boundaries. One must take some care in dealing with boundary traces. A distinguishing feature of the DGFE method is that the edges/faces of the partition Th play an important role in the formulation and analysis of the schemes. We define E I := e e = ∂K + ∩ ∂K − , ∃ K + , K − ∈ Th , measd−1 (∂K + ∩ ∂K − ) > 0 , where, for accurate bookkeeping, K + is the element with the larger index. It follows from conditions (T1) and (T2) on Th that C5 hK ± ≤ he ≤ hK ± , for some constant 0 < C5 ≤ 1 that is independent of h, where he := diamd−1 (e). In this case, it is common to write he ≈ hK + and he ≈ hK − . Suppose v ∈ H 2 (Th ), so that the zeroth and first boundary traces are welldefined on each edge (face) of each cell. We define the jump of v across e via [v] |e = v + |e − v − |e , where v + and v − denote the restrictions of v to K + and K − , respectively. Likewise, [∂n v] |e := [∇v] |e · n+ , where n+ is the outward-pointing unit normal of the element with the larger index, i.e., K + . For e ∈ E I we define the 1 + − average of v on e to be {v}|e := 2 v |e +v |e , and similarly, {∂n v}|e := {∇v}|e ·n+ . For any K ∈ Th and integer q ≥ 1, let Pq (K) denote the set of all polynomials of degree at most q on K. The discontinuous finite element space is defined via Y S (Th ) = Pq (Th ) := Pq (K) = v ∈ L2 (Ω) v|K ∈ Pq (K) . K∈Th

Clearly, S (Th ) ⊂ H 2 (Th ) ⊂ L2 (Ω). The elliptic DG bilinear form, αh (·, ·) on H 2 (Th ) × H 2 (Th ), corresponding to pure homogeneous Neumann boundary conditions, is defined as X αh (u, v) := (∇u, ∇v)K K∈Th

−

X

({∂n u}, [v])e + ([u], {∂n v})e − γh−1 ([u], [v]) e e ,

(12)

e∈E I

where γ is a positive constant (called the penalty parameter) that is chosen accordingly, and (·, ·)e is the d − 1 dimensional L2 inner product on e ∈ E I . For all


2217

v ∈ H 2 (Th ) define X

|||v|||2 :=

2

k∇vkK +

K∈Th

where kf ke =

p

X γ he 2 2 2 k[v]ke + k{∂n v}ke , he γ I

(13)

e∈E

(f, f )e , for any f ∈ L2 (e).

Lemma 2.1. The following basic facts are widely used in the DG literature: ˚2 (Th ) := H 2 (Th ) ∩ L2 (Ω). 1. ||| · ||| is a semi-norm on H 2 (Th ) and a norm on H 0 2. Continuity: There exists a positive constant Cα,1 such that, for all u, v ∈ H 2 (Th ), |αh (u, v)| ≤ Cα,1 |||u||| |||v|||. 3. Non-negativity: There exist positive constants γ0 and Cα,2 = Cα,2 (γ0 ), both independent of h, such that, for all γ ≥ γ0 , there holds αh (u, u) ≥ Cα,2 |||u|||2 , for all u ∈ S (Th ). 2 2 (Ω), 4. Consistency: If u ∈ HN p then − (∆u, v) = αh (u, v), for all v ∈ H (Th ). 5. Equivalency: |||v|||α := αh (v, v) is a semi-norm on S (Th ) and a norm on ˚ (Th ) := S (Th ) ∩ L2 . Moreover, for any v ∈ S (Th ) S 0 p p Cα,2 |||v||| ≤ |||v|||α ≤ Cα,1 |||v|||. (14) 6. Inverse Equivalency: If Th is globally quasi-uniform (T3), then for every v ∈ S (Th ), |||v|||α ≤ Ch−1 kvkL2 , (15) for some C > 0 that is independent of h. We needR a few technical results, which we now state. We use the notation u := |Ω|−1 Ω u dx, where |Ω| = measd (Ω), throughout the paper. The following basic result is well known. See [23, 24]. Lemma 2.2. Suppose the dimension of space is two (d=2) and Th is globally quasiuniform (T3). For any u ∈ H 2 (Th ) and any 2 ≤ p < ∞, and

ku − ukLp ≤ C|||u|||,

(16)

kukLp ≤ C |||u||| + kukL2 .

(17)

˚ (Th ), which is defined Later we will need the discrete Laplacian, ∆h : S (Th ) → S ˚ (Th ) denotes the unique solution to the as follows: for any vh ∈ S (Th ), ∆h vh ∈ S problem (∆h vh , χ) = −αh (vh , χ) , ∀ χ ∈ S (Th ) . (18) 2

In particular, setting χ = ∆h vh in (18), we obtain k∆h vh kL2 = −αh (vh , ∆h vh ). Proofs of the following two estimates involving the discrete laplacian can be found in [24]: Lemma 2.3. Suppose the dimension of space is two (d=2) and Th is globally quasiuniform (T3). Then, for any u ∈ S (Th ), 1

ku − ukL∞ k∇h ukL3

1

≤ C kukL2 2 k∆h ukL2 2 , 1 3

(19) 2 3

1 3

2 3

≤ C ku − ukL2 k∆h ukL2 + C kukL2 k∆h ukL2 .

(20)

2218


3. A mixed DG convex splitting scheme. 3.1. Definition of the scheme. We will assume from this point forward that θ ≥ 0. The sign of θ is a critical importance for the design of the scheme, as we explain below. Let M be a positive integer and 0 = t0 < t1 < · · · < tM = T be a uniform partition of [0, T ], with τ = ti −ti−1 , i = 1, . . . , M . Our mixed discontinuous Galerkin convex splitting scheme is defined as follows: for any 1 ≤ m ≤ M , given m um−1 ∈ S (Th ) find um h , wh ∈ S (Th ) such that h m m ∀χ ∈ S (Th ) , (δτ um h , χ) + αh (wh , χ) + θ (uh − u0 , χ) = 0, 1 3 m−1 m (um , ϕ + εαh (um ∀ϕ ∈ S (Th ) , h ) − uh h , ϕ) − (wh , ϕ) = 0, ε where m−1 um h − uh , u0h := Ph u0 , δτ um := h τ and Ph : H 2 (Th ) → S (Th ) is the elliptic projection:

(21a) (21b)

αh (Ph u − u, χ) = 0,

∀χ ∈ S (Th ) , (Ph u − u, 1) = 0. It is important to note that u0h − u0 , 1 = 0. We will make the very reasonable assumption that |u0 | ≤ 1. The elliptic projection is used in the initialization, mainly for simplicity in the forthcoming analysis, though other projections are possible. The scheme (21a) – (21b) is of convex splitting type and is expected to be of first order in time. There are two important properties that are automatically inherited in the convex splitting framework, unconditional energy stability, and unconditional unique solvability [16, 18, 28, 29, 32, 33]. We will demonstrate these properties later. Now, we can define a scheme that is equivalent to that above. For any 1 ≤ m ≤ ˚ (Th ), find v m ∈ S ˚ (Th ), wm ∈ S (Th ) such that M , given vhm−1 ∈ S h h (δτ vhm , χ) + αh (whm , χ) + θ (vhm , χ) = 0, ∀χ ∈ S (Th ) , 1 3 (vhm + u0 ) − vhm−1 + u0 , ϕ + εαh (vhm , ϕ) − (whm , ϕ) = 0, ε ∀ϕ ∈ S (Th ) ,

(22a)

(22b)

where vh0 := Ph u0 − u0 .

Hence, vh0 , 1 = 0. By setting χ ≡ 1 in (21a) and (22a) and observing that αh (1, ϕ) = 0 for all ϕ ∈ S (Th ), one finds that, provided solutions for the two schemes exist, they are related via vhm + u0 = um h ,

˚ (Th ) , vhm ∈ S

˚ um h ∈ S (Th ) + u0 ,

for all 1 ≤ m ≤ M . In either case, it is generally true that the average mass of whm will change with the time step m, i.e., (whm , 1) 6= whm−1 , 1 , in general. 3.2. Unconditional solvability. In this subsection we show that our fully discrete scheme is unconditionally uniquely solvable. We begin by defining an invertible ˚ (Th ) → S ˚ (Th ) via the following variational problem: given linear operator Th : S ˚ ˚ ζ ∈ S (Th ), find Th (ζ) ∈ S (Th ) such that αh (Th (ζ), χ) = (ζ, χ)

˚ (Th ) . ∀χ ∈ S

(23)


2219

˚ (Th ). This clearly has a unique solution because αh ( · , · ) is an inner product on S We now wish to define a mesh-dependent “−1” norm, i.e., a discrete analogue to the H −1 norm. We omit the details of the next result, the discrete analog of Lemma 1.1, for brevity. ˚ (Th ) and set Lemma 3.1. Let ζ, ξ ∈ S (ζ, ξ)−1,h := αh (Th (ζ), Th (ξ)) = (ζ, Th (ξ)) = (Th (ζ), ξ) .

(24)

˚ (Th ), and the induced negative norm is ( · , · )−1,h defines an inner product on S equivalent to the operator norm: q (ζ, χ) sup kζk−1,h := (ζ, ζ)−1,h = . (25) |||χ||| α ˚ h) 06=χ∈S(T ˚ (Th ), Consequently, for all χ ∈ S (Th ) and all ζ ∈ S |(ζ, χ)| ≤ kζk−1,h |||χ|||α .

(26)

The following Poincare-type estimate holds: kζk−1,h ≤ C kζkL2 ,

˚ (Th ) , ∀ζ∈S

(27)

for some C > 0 that is independent of h. Finally, if Th is globally quasi-uniform (T3), then the following inverse estimate holds: kζkL2 ≤ Ch−1 kζk−1,h ,

˚ (Th ) , ∀ζ∈S

(28)

for some C > 0 that is independent of h. Using our discrete negative norm we can define a variational problem closely related to our fully discrete scheme. ˚ (Th ) be given. Set β := 1 + τ θ. For all vh ∈ S ˚ (Th ), Lemma 3.2. Let vhm−1 ∈ S define the nonlinear functional

2 βvh − vhm−1 τ 1 4

Gh (vh ) := + kvh + u0 kL4

2β τ 4ε −1,h ε 1 + |||vh |||2α − vhm−1 + u0 , vh . 2 ε ˚ (Th ). Consequently, Gh Gh is strictly convex and coercive on the linear subspace S m ˚ ˚ (Th ) is the unique has a unique minimizer, call it vh ∈ S (Th ). Moreover, vhm ∈ S minimizer of Gh if and only if it is the unique solution to 1 m−1 1 m 3 m ,ϕ = (29) (vh + u0 ) , ϕ + εαh (vhm , ϕ) − wh,? v + u0 , ϕ ε ε h ˚ (Th ), where wm ∈ S ˚ (Th ) is the unique solution to for all ϕ ∈ S h,? m vh − vhm−1 m ˚ (Th ) . αh wh,? , χ = − , χ − θ (vhm , χ) ∀χ ∈ S (30) τ Proof. The strict convexity of Gh is easily established as in [18], and we omit this step. The coercivity of Gh follows from an easy-to-establish estimate of the form ˚ (Th ) , Gh (uh ) ≥ M1 |||uh |||2α − M2 ∀uh ∈ S where 0 < M1 , M2 < ∞ are constants that depend on ε and θ. By the standard ˚ (Th ), theory of convex optimization [11], Gh has a unique (bounded) minimizer in S

2220


call it vhm . Moreover, vhm is the unique minimizer of Gh if and only if it is the unique ˚ (Th ), where solution to the Euler-Lagrange equation hA (vhm ) , ϕi = 0, for all ϕ ∈ S d hA (vh ) , ϕi := Gh (vh + sϕ) ds s=0 1 3 = (vh + u0 ) , ϕ + εαh (vh , ϕ) ε βvh − vhm−1 1 m−1 ,ϕ − + vh + u0 , ϕ . τ ε −1,h To conclude, we define m wh,? := −Th

βvhm − vhm−1 τ

˚ (Th ) . ∈S

Finally, we are in the position to prove the unconditional unique solvability of our scheme. Theorem 3.3. The scheme (21a) – (21b), or, equivalently, the scheme (22a) – (22b), is uniquely solvable for any mesh parameters τ and h and for any phase parameter ε. Proof. Suppose vhm−1 , 1 = 0. It is clear that a necessary condition for solvability of (22a) – (22b) is that (vhm , 1) = 0, as can be found by taking χ ≡ 1 in (22a). Now, m ˚ (Th ) × S ˚ (Th ) be a solution of (29) – (30). Set let vhm , wh,? ∈S wm h :=

u0 1 1 (vhm + u0 )3 − (vhm + u0 ) , 1 = (vhm + u0 )3 , 1 − , ε|Ω| ε|Ω| ε

m and define whm := wh,? +wm h . There is a one-to-one correspondence of the respective m m ˚ ˚ (Th ) is a solution to (29) – (30) if and only if solution sets: vh , wh,? ∈ S (Th ) × S ˚ (Th ) × S (Th ) is a solution to (22a) – (22b) if and only if um , wm ∈ vhm , whm ∈ S h h S (Th ) × S (Th ) is a solution to (21a) – (21b), where m um h = vh + u0 ,

m whm = wh,? + wm h .

But (29) – (30) admits a unique solution, which proves that (21a) – (21b) and (22a) – (22b) are uniquely solvable. 3.3. Unconditional energy stability. We now show that the solutions to our scheme enjoy stability properties that are similar to those of the PDE solutions, and, moreover, these properties hold regardless of the sizes of h and τ . The first property, the unconditional energy stability, is a direct result of the convex decomposition represented in the scheme. m Lemma 3.4. Let (um h , wh ) ∈ S (Th ) × S (Th ) denote the unique solution of the scheme (21a) – (21b). Then the following energy law holds for any h, τ > 0: ( ` ` X X

1 ε 2 2 2 2 2

δτ (um |||δτ um Eh u`h + τ kδτ um h k−1,h + τ h |||α + h ) L2 2 4ε m=1 m=1 (31) ) 1 θ 1 m m 2 m 2 m 2 0 + ku δτ uh kL2 + kδτ uh kL2 + kδτ uh k−1,h = Eh uh , 2ε h 2ε 2


for all 0 ≤ ` ≤ M , where 1 Eh (uh ) := 4ε

2221

the spatially discrete energy is given by

2

ε θ

2 2

(uh ) − 1 2 + |||uh |||2α + kuh − u0 k−1,h . 2 2 L

(32)

m Proof. We first set χ = Th (δτ um h ) in (21a) and set ϕ = δτ uh in (21b) to obtain 2

m m m m kδτ um = 0, h k−1,h + (wh , δτ uh ) + θ (uh − u0 , δτ uh )−1,h 1 3 m−1 m m m (um + εαh (um , δτ um h ) − uh h , δτ uh ) − (wh , δτ uh ) = 0. h ε Adding (33) and (34), using the identities 1 m 2 m 2 αh (um δτ |||um h , δτ uh ) = h |||α + τ |||δτ uh |||α , 2

2 1

τ h

2 m−1 2 2 m 3 m

δτ (um (uh ) − uh , δτ uh = δτ (um h ) − 1 2 + h ) L2 4 4 L i 2

m (um h − u0 , δτ uh )−1,h

and applying the operator τ

(33) (34)

2

m m + 2 kum h δτ uh kL2 + 2 kδτ uh kL2 , i 1h 2 m 2 = δτ kum u k + τ kδ u k − 0 τ h −1,h , h −1,h 2

P`

m=1

to the combined equation the result is obtained.

The discrete energy law immediately implies the following uniform (in h and τ ) m a priori estimates for um h and δτ uh . m Corollary 3.5. Let (um h ,w h ) ∈ S (Th ) × S (Th ) be the unique solution of (21a)– 0 (21b). Suppose that E uh < C, independent of h. Then the following estimates hold for any h, τ > 0:

2

m 2

2 m 2 m 4 m 2 + |||u ||| + ) − 1 + ku k max kum − u k + ku k

(uh

0 −1,h h α ≤ C, h L4 h h L2 L2

0≤m≤M

(35) τ

M X

2

kδτ um h k−1,h ≤ C,

m=1

(36) for some constant C > 0 that is independent of h, τ , and T . The next set of a priori stability estimates are completely unconditional; they are proven without any restrictions on h, τ , or ε. m Lemma 3.6. Suppose d = 2 and Th is globally quasi-uniform (T3). Let (um h ,wh ) ∈ 0 S (Th ) × S (Th ) be the unique solution of (21a)–(21b). Suppose that E uh < C,

0 2

w 2 < C independent of h, where w0 is defined below in (42), and that T ≥ 1 h L

h

(for simplicity). The following estimates hold for any h, τ > 0: M X 2 m 2 m 2 m 4 + kw k + kδ u k + ku k ≤ T C, τ |||whm |||2α + k∆h um k τ h L2 h L2 h L2 h L∞

(37)

m=1

max

1≤m≤M

2

2

4

m kwhm kL2 + k∆h um h kL2 + kuh kL∞

for some constant C > 0 that is independent of h, τ , and T .

≤ T C,

(38)

2222


Proof. Setting χ = whm in (21a) and using (26), (35) and (36), we have m m m |||whm |||2α = αh (whm , whm ) = − (δτ um h , wh ) − θ (uh − u0 , wh )

1 2 2 m m 2 ≤C kδτ um h k−1,h + C kuh − u0 k−1,h + |||wh |||α 2 1 2 m 2 ≤C kδτ um h k−1,h + |||wh |||α + C. 2

(39)

PM Applying τ m=1 gives (37.1) – which, in our notation, is the bound on the first term of the left side of (37). m Setting ϕh = ∆h um h in (21b) and using the definition of ∆h uh , we get 2

1 m 3 (uh ) − um−1 , ∆h um h h ε

2 1 1 ε

m 3 2 m−1 m 2 2 ≤ |||wh |||α + |||um k∆h um

2. h |||α + h kL2 + C (uh ) − uh 2 2 2 L

m m ε k∆h um h kL2 = − (wh , ∆h uh ) +

Hence,

2

m 3 2 m−1 m 2 m 2 ε k∆h um k ≤ |||w ||| + |||u ||| + C ) − u

(u

2. 2 h L h α h α h h

(40)

L

Now using (35), and (17), we have

m−1 2

m 3

6

u

2 ≤ C.

(uh ) − um−1 2 2 ≤ 2 kum + k 6 h L h h L L

(41)

Putting the last two inequalities together, we have 2

m 2 ε k∆h um h kL2 ≤ |||wh |||α + C

PM Applying τ m=1 , estimate (37.2) follows. Now, take ϕ = whm in (21b). Then, using (35) and (41) and rearranging terms, we get 2

kwhm kL2 ≤ C + ε|||whm |||2α . PM Applying τ m=1 , estimate (37.3) now follows from (37.1). Next we prove (37.4) and (38.1) together. To do so, we first define wh0 via 1 0 3 wh0 , ϕ := εαh u0h , ϕ + uh − u0h , ϕ , ε

(42)

for all ϕ ∈ S (Th ). We also define δτ u0h :≡ 0. Now, we subtract (21b) from itself at consecutive time steps to obtain 1 m 3 m−1 3 (u ) − u , ϕ whm − whm−1 , ϕ =τ εαh (δτ um , ϕ) + h h h ε τ − δτ um−1 ,ϕ , h ε

(43)


2223

for all ϕ ∈ S (Th ), which is well-defined for all 1 ≤ m ≤ M . Taking ϕ = whm in (43) and χ = τ εδτ um h in (21a) and adding the results yields 2

τ (δτ whm , whm ) + τ ε kδτ um h kL2 m τ m−1 2 m 2 m m−1 = δτ um (u ) + u u + u , wh h h h h h ε τ m δτ um−1 − εθτ (um , whm h − u0 , δτ uh ) − h ε τ

2 m−1 2 m m−1 ) + u u + u ≤ (um

4 kwhm kL4 kδτ um h h h h kL2 h ε L

τ m−1 m m

+ εθτ |||um , h |||α kδτ uh k−1,h + |||wh |||α δτ uh −1,h ε

2 2

2 m m−1 m 2 m 2 ≤τ C (um + um−1 kw k + |||w |||

2 h ) + uh uh h h α L h L4

2 τε m−1 m 2 m 2 kδτ uh kL2 + τ C 1 + |||wh |||α + δτ uh −1,h , + 2 where we have used (17), (26) from Lemma 3.1, and (35). To continue, we use (35) and (16) to obtain 2

τ (δτ whm , whm ) + τ ε kδτ um h kL2

m−1 4 m 2 4

8 kwh k 2 + |||whm |||2α ≤Cτ kum h kL8 + uh L L

2 τε m 2 m 2 + kδτ uh kL2 + τ C 1 + |||wh |||α + δτ uhm−1 −1,h 2 2 m−1 4 4 ≤Cτ u40 + |||um |||α kwhm kL2 + |||whm |||2α h |||α + |||uh

τε 2 m−1 2 m 2

+ kδτ um h kL2 + τ C 1 + |||wh |||α + δτ uh −1,h τε 2 2 2 ≤Cτ kwhm kL2 + |||whm |||2α + kδτ um h kL2 2

2 2

. + τ C 1 + |||whm |||2α + δτ um−1 h −1,h Then, τ (δτ whm , whm ) +

τε 2 m−1 2 m 2 m 2

kδτ um . h kL2 ≤ Cτ kwh kL2 + |||wh |||α + Cτ δτ uh −1,h 2

Summing from m = 1 to m = `, and using (36), (37.1), (37.3), δτ u0h ≡ 0, we have ` X

1 2

wh` 2 2 − 1 wh0 2 2 + τ ε kδτ um h kL2 L L 2 2 2 m=1

≤Cτ

` X

` X

2

δτ um−1 2 kwhm kL2 + |||whm |||2α + Cτ + CT ≤ CT. h −1,h

m=1

m=1

Finally, for any 1 ≤ ` ≤ M , ` X

1 1 2

wh` 2 2 + τ ε

wh0 2 2 ≤ CT, kδτ um h kL2 ≤ CT + L L 2 2 m=1 2

2 where we used the assumption that wh0 L2 is bounded above independent of h in the last step. This proves (37.4) and (38.1).

2224


m Again, setting ϕh = ∆h um h in (21b) and using the definition of ∆h uh , 2

m m ε k∆h um h kL2 = −εαh (uh , ∆h uh )

1 m 3 m (uh ) − um−1 , ∆ u h h h ε

2 ε ε 1 m 3 2 2 m−1 ) − u + k∆h um k +

(u

2 + k∆h um h h L2 h kL2 . h 4 ε2 4 L

= − (whm , ∆h um h )+ ≤

1 2 kwhm kL2 ε

Consequently, using (37.3) and (41), 2

ε k∆h um h kL2 ≤

2 2 2

3 2 m−1 ) − u kwhm kL2 + 2 (um

2 ≤ T C, h h ε ε L

giving estimate (38.2). To prove estimate (37.5), consider 1

1

m m 2 m 2 kum h kL∞ = kvh + u0 kL∞ ≤ |u0 | + C kvh kL2 k∆h vh kL2 ,

where we have used the globally quasi-uniformity of Th and (19). This leads to 4

2

2

2

m m m kum h kL∞ ≤ C + C kvh kL2 k∆h vh kL2 ≤ C + C k∆h uh kL2 .

(44)

PM Applying τ m=1 and using (37.2), estimate (37.5) follows. Estimate (38.3), follows from (44) and (38.2). Remark 3.7. A version of the last set of a priori stability estimates is given in [24] for a different scheme applied to a similar model equation. However, unlike in [24], we are able to prove these stability estimates without any restrictions of h and τ , and, moreover, the bounds in our estimates grow only linearly in T . These points represent important advantages of our approach.

4. Error estimates for the fully discrete convex splitting scheme. For the error estimates that we pursue in this section, we assume that Ω ⊂ R2 is convex polygonal and that weak solutions (u, w) to (2a) − (2b) have the regularities 2 u0 ∈ H q+1 (Ω) ∩ HN (Ω), 1,∞ q+1 2,∞ 2 1 2 u∈W (0, T ; H (Ω)) ∩ W (0, T ; H (Ω)) ∩ H 0, T ; H (Ω) , w ∈ L∞ (0, T ; H q+1 (Ω)) ∩ L2 0, T ; H 2 (Ω) .

(45)

In this case, (u, w) solve the following broken variational problem: (∂t u, χ) + αh (w, χ) + θ (u − u0 , χ) = 0, 1 3 u − u, ϕ + εαh (u, ϕ) − (w, ϕ) = 0, ε

∀χ ∈ H 2 (Th ) ,

(46a)

∀ϕ ∈ H 2 (Th ) .

(46b)

Define the time lag operator, Lτ u(t) := u(t−τ ), the backward difference operator, δτ u(t) := τ −1 (u(t) − Lτ u(t)), and the following error terms: for f = u, w and 0 ≤ t ≤ T, f EA (t) := f (t) − Ph f (t),

σhu := δτ Ph u(t) − ∂t Ph u(t).

(47)


2225

Then, from (46a) and (46b) we have the following: for all χ, φ ∈ S (Th ), u (δτ Ph u, χ) + αh (Ph w, χ) + θ (Ph u − u0 , χ) = (σhu , χ) − (∂t EA , χ) u − θ (EA , χ) , (48a) τ w , ϕ) εαh (Ph u, ϕ) − (Ph w, ϕ) = (δτ u, ϕ) + (EA ε 1 3 u − Lτ u, ϕ . (48b) − ε Now, we define the piecewise constant approximations: for m = 1, . . . M and for t ∈ (tm−1 , tm ], u ˆ( · , t) := um w( ˆ · , t) := whm ( · ), h ( · ), m m where uh and wh are the solutions of the fully discrete convex splitting scheme (21a) – (21b). Thus, for all χ, φ ∈ S (Th ),

(δτ u ˆ, χ) + αh (w, ˆ χ) + θ (ˆ u − u0 , χ) = 0, εαh (ˆ u, ϕ) − (w, ˆ ϕ) = − For f = u, w and 0 ≤ t ≤ T , setting E f (t) := Ph f (t) − fˆ(t) ∈ S (Th ) , h

(49a) 1 3 u ˆ − Lτ u ˆ, ϕ . ε

(49b)

E f (t) := f (t) − fˆ(t) ∈ H 2 (Th ) ,

and subtracting (49a) from (48a) and (49b) from (48b), we have, for all χ, φ ∈ S (Th ), u u (δτ Ehu , χ) + αh (Ehw , χ) + θ (Ehu , χ) = (σhu , χ) − (∂t EA , χ) − θ (EA , χ) , (50a) τ w εαh (Ehu , ϕ) − (Ehw , ϕ) = (δτ u, ϕ) + (EA , ϕ) ε 1 3 − u −u ˆ3 − Lτ E u , ϕ . (50b) ε Setting χ = Th (δτ Ehu ) in (50a) and ϕ = δτ Ehu in (50b) and by adding the resulting equations we obtain 2

εαh (Ehu , δτ Ehu ) + kδτ Ehu k−1,h +θ (Ehu , δτ Ehu )−1,h u u = (σhu − ∂t EA − θEA , Th (δτ Ehu )) τ (51) w + (δτ u, δτ Ehu ) + (EA , δτ Ehu ) ε 1 3 − u −u ˆ3 − Lτ E u , δτ Ehu . ε The last expression is the key error equation. We now proceed to estimate the terms on the right hand side of (51).

Lemma 4.1. Suppose that (u, w) is a weak solution to (46a) – (46b) with the regularities (45). Then, for any h, τ > 0, there exists C > 0, independent of h and τ , such that u u 2 kσhu − ∂t EA − θEA kL2 ≤ C h2q+2 + τ 2 . (52) Proof. Taylor’s theorem and well-known properties of the elliptic projection imply the estimates kσhu ( · , t)kL2 ≤ Cτ,

u k∂t EA ( · , t)kL2 ≤ Chq+1 ,

u kEA ( · , t)kL2 ≤ Chq+1 ,

for t ∈ [0, T ], and the result follows. We need a technical lemma that will be used a number of times.

2226


Lemma 4.2. Suppose Th is globally quasi-uniform (T3), g ∈ H 2 (Th ), and v ∈ ˚ (Th ). Then S |(g, v)| ≤ C|||g||| kvk−1,h , (53) for some C > 0 that is independent of h. Proof. If g ∈ S (Th ), we can apply Lemma 3.1 directly and use the equivalence of ||| · ||| and ||| · |||α as semi-norms on S (Th ). Otherwise, using the triangle inequality, the Cauchy-Schwarz inequality, and Lemma 3.1, |(g, v)| ≤ |(g − Ph g, v)| + |(Ph g, v)| ≤ kg − Ph gkL2 kvkL2 + |||Ph g|||α kvk−1,h . Using the standard elliptic projection estimate kg − Ph gkL2 ≤ Ch|||g|||, we have |(g, v)| ≤ Ch|||g||| kvkL2 + |||Ph g|||α kvk−1,h . Finally, using the (uniform) inverse estimate h kvkL2 ≤ C kvk−1,h from Lemma 3.1, and the stability of the elliptic projection, |||Ph g|||α ≤ C|||g|||, we have the result. Lemma 4.3. Suppose that d = 2 and Th is globally quasi-uniform (T3). Let (u, w) be a weak solution to (46a) – (46b) with the regularities (45). Then, for any h, τ > 0, |||u3 ( · , t) − u ˆ3 ( · , t)||| ≤ C|||E u ( · , t)|||, (54) for a.e. t ∈ [0, T ], where E u := u − u ˆ. Proof. Since the estimates are tedious, we suppress the details for brevity. One begins by writing u3 − u ˆ3 = E u u2 + uˆ u+u ˆ2 . (55) Then, using the definition of the norm ||| · ||| in (13), the unconditional a priori estimates in Lemma 3.6, specifically (38.2) and (38.3), the estimates (16) and (20), and standard trace inequalities, the result follows. Remark 4.4. A version of this result is proven in [24, Lemma 3.2] under the simplifying (but not limiting) assumption that u and its DG approximation are mean-zero functions. We point out that, because their a priori energy estimates – the counterparts of (38.2) and (38.3) in [24] – come with time step restrictions, Lemma 3.2 in [24] does not hold unconditionally, in contrast with the result here. Remark 4.5. From here we assume that d = 2 and Th is globally quasi-uniform (T3) in order to make use of Lemmas 4.2 and 4.3. Lemma 4.6. Suppose that (u, w) is a weak solution to (46a) – (46b) with the regularities (45). Then, for any h, τ > 0, there exists a constant C > 0, independent of h and τ , such that εαh (Ehu , δτ Ehu ) +

1 2 kδτ Ehu k−1,h + θ (Ehu , δτ Ehu )−1,h 2 ≤ C|||Ehu |||2α + C|||Lτ Ehu |||2α + C(τ 2 + h2q ).

(56)


2227

Proof. Using Lemma 4.1, the Cauchy-Schwarz inequality, and the estimate 2

2

kTh (δτ Ehu )kL2 ≤ C|||Th (δτ Ehu ) |||2α = C kδτ Ehu k−1,h , we get the following estimate immediately: u u |(σhu − ∂t EA − θEA , Th (δτ Ehu ))| ≤ C(τ 2 + h2q+2 ) +

1 2 kδτ Ehu k−1,h . 10

(57)

Now, from the standard DG approximation theory w |||EA ||| = |||Ph w − w||| ≤ Chq .

Applying Lemma 4.2 and the last estimate w w |(EA , δτ Ehu )| ≤C|||EA ||| kδτ Ehu k−1,h

≤Chq kδτ Ehu k−1,h (58) 1 u 2 2q kδτ Eh k−1,h . ≤Ch + 10 Now observe that, since u ∈ W 2,∞ 0, T ; H 2 (Ω) , it follows that τ |||δτ u||| ≤ Cτ . Consequently, using Lemma 4.2, we have τ |(δτ u, δτ Ehu )| ≤Cτ |||δτ u||| kδτ Ehu k−1,h ε (59) 1 2 ≤Cτ 2 + kδτ Ehu k−1,h . 10 u u Using Lemmas 4.2 and 4.3, as well as E = EA + Ehu and a standard DG error estimate, 1 3 u −u ˆ3 , δτ Ehu ≤C|||u3 − u ˆ3 ||| kδτ Ehu k−1,h ε 1 2 ≤C|||u3 − u ˆ3 |||2 + kδτ Ehu k−1,h 10 1 2 (60) ≤C|||E u |||2 + kδτ Ehu k−1,h 10 1 2 u 2 ≤C|||EA ||| + C|||Ehu |||2α + kδτ Ehu k−1,h 10 1 2 ≤Ch2q + C|||Ehu |||2α + kδτ Ehu k−1,h . 10 Finally, with similar steps as in the last estimate, 1 |(Lτ E u , δτ Ehu )| ≤C|||Lτ E u ||| kδτ Ehu k−1,h ε ≤Ch2q + C|||Lτ Ehu |||2α +

1 2 kδτ Ehu k−1,h . 10

(61)

Combining estimates (57) – (61) with the error equation (51) and using the triangle inequality, the result follows. Lemma 4.7. Suppose that (u, w) is a weak solution to (46a) – (46b) with the regularities (45). Then, for any h, τ > 0, there exists a constant C > 0, independent of h and τ , such that 2 2 |||Ehw |||2α ≤ C kδτ Ehu k−1,h + C kEhu k−1,h + C h2q+2 + τ 2 . (62)

2228


Proof. Define Ehw : [0, T ] → R via 1 1 Ehw (t) := (E w (t), 1) = (Ph w(t) − w(t), ˆ 1) . |Ω| h |Ω|

(63)

u u Since (σhu (t) − ∂t EA (t) − θEA (t), 1) = 0, for a.e. t ∈ [0, T ], setting χ = Ehw in (50a), we have u u |||Ehw |||2α = − (δτ Ehu + θEhu , Ehw ) + (σhu − ∂t EA − θEA , Ehw ) u u = − (δτ Ehu + θEhu , Ehw ) + σhu − ∂t EA − θEA , Ehw − Ehw (t)

u u ≤ kδτ Ehu + θEhu k |||Ehw |||α + kσhu − ∂t EA − θEA k 2 Ehw − E w (t) 2 −1,h

L

h

L

u u ≤ kδτ Ehu + θEhu k−1,h |||Ehw |||α + C kσhu − ∂t EA − θEA kL2 |||Ehw |||α 1 2 2 ≤C kδτ Ehu k−1,h + C kEhu k−1,h + C h2q+2 + τ 2 + |||Ehw |||2α , 2 where we have used (16) and (26). The result now follows.

Using the last two lemmas, we are ready to show the main convergence result for our convex splitting scheme. Theorem 4.8. Suppose that (u, w) is a weak solution to (46a) – (46b) with the regularities (45). Then, provided 0 < τ < τ0 , for some τ0 sufficiently small, M h h i i X 2 2 max |||Ehu (tm )|||2α + kEhu (tm )k−1,h + τ |||Ehw (tm )|||2α + kδτ Ehu (tm )k−1,h 1≤m≤M

m=1

≤ C(T )(τ 2 + h2q ), (64) for some C(T ) > 0 that is independent of τ and h. Proof. Setting t = tm and using Lemma 4.6 and the arithmetic-geometric mean inequality, we have ε 1 θ 2 2 δτ |||Ehu (tm )|||2α + kδτ Ehu (tm )k−1,h + δτ kEhu (tm )k−1,h 2 2 2 ≤ C|||Ehu (tm )|||2α + C|||Lτ Ehu (tm )|||2α + C(τ 2 + h2q ). P` Let 1 < ` ≤ M . Applying τ m=1 and using Ehu (t0 ) ≡ 0, ` θ ε τ X 2 2 |||Ehu (t` )|||2α + kδτ Ehu (tm )k−1,h + kEhu (t` )k−1,h 2 2 m=1 2 2

` X

2q

≤ CT (τ + h ) + C1 τ

(65)

|||Ehu (tm )|||2α .

m=1

If 0 < τ ≤ τ0 :=

ε 4C1

0, the energy term θ2 ku − u0 kH −1 in (8) is convex and stabilizing and its presence in the energy may lead to some suppression of the coarsening process during phase decomposition. This can be seen in Fig. 1, columns 1 and 2. Another interesting feature of this case is the possible presence of stable isolated precipitate islands that persist throughout phase decomposition and coarsening [9, 10]. A complete arrest of the coarsening process requires that θ is sufficiently large. For the DGFE simulation in Fig. 5, we use θ = 15000.0 and we take the average composition u0 = −0.3 (which is rather far from the symmetric point u0 = 0.0). The other parameters are as follows: √ Ω = (0, 8)2 , ε = 0.02, h = 81282 , and τ = 2.0 × 10−5 . The mesh is like that shown in Fig. 2, but with seven levels of global quadrisection refinement used to obtain the final uniform mesh. For the simulation in Fig. 5a, initially there are two isolated regions – one near x = 2, y = 0, and one near x = 6, y = 0 – where small amplitude perturbations about the average u0 are introduced. Elsewhere u = u0 exactly. The perturbations rapidly grow in amplitude, then propagate radially outward. The result is a pseudocrystallization of the domain, similar to that predicted by the phase field crystal (PFC) equation [6, 22]. Note that, in this case, θ is sufficiently large to prevent the coarsening of the atom-like precipitate particles and a crystalline lattice of islands will persist over long times. This is in contrast to the typical situation modeled by the Cahn-Hilliard equation (θ = 0.0), where coarsening dominates long time


2235

dynamics. The discrete energy for this simulation is plotted in Fig. 5(b). We observe that, as predicted in Lemma 3.4, the discrete energy is non-increasing in time. t = 0.0

t = 0.0004

t = 0.0008 140 120

t = 0.0066

t = 0.015

t = 0.02

Energy

100 80 60 40 20 0

(a)

0

0.005

0.01 Time

0.015

0.02

(b)

Figure 4. (Color figure.) (a) Snapshots from a simulation of spinodal decomposition using quadratic elements, q = 2. The domain Ω is a dodecagon of diameter 3.2, and we use a uniform triangulation consisting of 12288 elements. The base level mesh consists of 12 triangles, and there are exactly 5 levels of global quadrisection refinement to obtain the final mesh. The final time is T = 0.02, and τ = 10−6 . The interface parameter is ε = 0.0125; u0 ≈ 3.0 × 10−03 ; and θ = 0.0. (b) The discrete Cahn-Hilliard energy (32) corresponding to the simulation in (a), plotted as a function of time. We observe that, as predicted in Lemma 3.4, the energy is nonincreasing in time.

6.3. Nonlinear multigrid solver convergence test. We now show some evidence suggesting that our multigrid solver, described in Section 5, has near optimal complexity. Residuals (precisely, the 2-norms of the algebraic residual vectors) for our FAS nonlinear multigrid solver are plotted against the number of V-cycle iterations in Fig. 6. For the tests we use quadratic, q = 2, elements with the same problem as in the convergence test case where (71) is the exact solution. The parameters are Ω = (0, 1)2 ; τ = 1.5 × 10−1 ; θ = 5.0, and ε = 0.5. We used `max = 2, i.e., 2 pre- and 2 post-smoothing passes in multigrid parlance [27]. We plot the residuals only at time t = 1.5, i.e., after exactly 10 time steps, and vary the space √ √ step sizes from h = 22 to h = 2562 , in multiples of 12 , but keep τ fixed. The trian√ √ √ √ √ gulations corresponding to h = 2, h = 22 , h = 42 , h = 82 , and h = 162 are as shown in Fig. 2. We observe that the residual decreases by nearly the same factor for each V-cycle iteration, thus showing that the solver is essentially independent of h. More complete analysis and performance testing for the solver are planned for a latter paper. Acknowledgments. Some of this work was completed after AA took up residence at the Statistical and Applied Mathematical Sciences Institute (SAMSI). AA gratefully acknowledge partial financial support from the National Science Foundation under Grant DMS-1127914 to that institute. OK and SMW gratefully acknowledge

2236

A. C. ARISTOTELOUS, O. KARAKASHIAN AND S. M. WISE t = 0.2 x 10-02

t = 0.6 x 10-02

664 662 660

t = 1.4 x 10-02

Energy

658 t = 1.0 x 10-02

656 654 652 650 648

0

0.002 0.004 0.006 0.008

(a)

0.01 Time

0.012 0.014 0.016 0.018

0.02

(b)

Figure 5. (Color figure.) (a) Snapshots from the simulation of phase decomposition using quadratic elements, q = 2. The parameters are as follows Ω = (0, 8)2 , ε = 0.02, θ = 15000.0, u0 = −0.3, √ 8 2 h = 128 , and τ = 2 × 10−5 . The mesh is like that shown in Fig. 2 but with seven levels of refinement. At t = 0 there are two regions – one near x = 2, y = 0, and one near x = 6, y = 0 – where small amplitude perturbations of the solution are introduced. The perturbations rapidly grow in amplitude, then propagate radially outward and the domain “crystalizes.” (b) The discrete energy (32), corresponding to the simulation in (a), plotted as a function of time. We observe that, as predicted in Lemma 3.4, the energy is non-increasing in time. The energy levels off noticeably after about t = 0.0155. This corresponds to the time when the entire domain is phase separated, i.e., when the dots reach the top of the computational domain. partial financial support from the National Science Foundation through grants DMS 0811314 and DMS 1115390, respectively. REFERENCES [1] A. Aristotelous, “Adaptive Discontinuous Galerkin Finite Element Methods for a Diffuse Interface Model of Biological Growth,” PhD thesis, University of Tennessee, 2011. [2] D. N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM Journal on Numerical Analysis, 19 (1982), 742–760. [3] D. N. Arnold, F. Brezzi, B. Cockburn and L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis, 39 (2002), 1749–1779. [4] I. Babuˇska and M. Zl´ amal, Nonconforming elements in the finite element method with penalty, SIAM Journal on Numerical Analysis, 10 (1973), 863–875. [5] G. A. Baker, Finite element methods for elliptic equations using nonconforming elements, Math. Comp., 31:45–59, 1977. [6] A. Baskaran, Z. Hu, J. S. Lowengrub, C. Wang, S. M. Wise and P. Zhou, Energy stable and efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal equation, Journal of Computational Physics, 250 (2013), 270–292. [7] J. H. Bramble, “Multigrid Methods,” Research Notes in Mathematics Series. Chapman and Hall/CRC, London, 1993. [8] J. W. Cahn, On spinodal decomposition, Acta Metallurgica, 9 (1961), 795–801.


0.01

h=d/2 h=d/4 h=d/8 h=d/16 h=d/32 h=d/64 h=d/128 h=d/256

0.0001 1e-06 log of residual

2237

1e-08 1e-10 1e-12 1e-14 1e-16

0

2

4

6 8 number of V cycles

10

12

14

Figure 6. Residuals for the FAS multigrid solver plotted against the number of V-cycle iterations. We use quadratic (q = 2) elements with the same problem as in the convergence test case where (71) is the exact solution. The parameters are Ω = (0, 1)2 ; τ = 0.15; θ = 5.0; and ε = 0.5. We plot the residuals only at time t = 1.5, i.e., after exactly 10 time steps. We vary the√space step size, as indicated, but keep τ fixed. In the figure, d = 2. We observe that the decrease in the residual after each V-cycle iteration is nearly independent of h. The hierarchical mesh used for the test is like that show in Fig. 2.

[9] R. Choksi, M. Maras and J. F. Williams, 2D phase diagram for minimizers of a Cahn-Hilliard functional with long-range interactions, SIAM Journal on Applied Dynamical Systems, 10 (2011), 1344–1362. [10] R. Choksi and X. Ren, On a derivation of a density functional theory for microphase separation of diblock copolymers, Journal of Statistical Physics, 113 (2003), 151–176. [11] P. G. Ciarlet, “Introduction to Numerical Linear Algebra and Optimisation,” Cambridge University Press, Cambridge, UK, 1989. [12] C. Collins, J. Shen and S. M. Wise, Unconditionally stable finite difference multigrid schemes for the Cahn-Hilliard-Brinkman equation, Commun. Comput. Phys., 13 (2013), 929–957. [13] J. Douglas and T. Dupont, Interior penalty procedures for elliptic and parabolic Galerkin methods, In “Computing Methods in Applied Sciences,” pages 207–216. Springer, Berlin, 1976. [14] C. M. Elliott and S. Zheng, On the Cahn-Hilliard equation, Archive for Rational Mechanics and Analysis, 96 (1986) ,339–357. [15] C. M. Elliott, The Cahn-Hilliard model for the kinetics of phase separation, In J.F. Rodrigues, editor, Mathematical Models for Phase Change Problems: Proceedings of the Euro´ pean Workshop held at Obidos, Portugal, October 1-3, 1988, International Series of Numerical Mathematics, 35–73, Berlin, 1989. Birkh¨ auser Verlag.

2238


[16] D. Eyre, Unconditionally gradient stable time marching the Cahn-Hilliard equation, In J. W. Bullard, R. Kalia, M. Stoneham, and L.Q. Chen, editors, Computational and Mathematical Models of Microstructural Evolution, volume 53, pages 1686–1712, Warrendale, PA, USA, 1998. Materials Research Society. [17] X. Feng and O. A. Karakashian, Fully discrete dynamic mesh discontinuous Galerkin methods for the Cahn-Hilliard equation of phase transition, Math. Comput., 76 (2007), 1093–1117. [18] X. Feng and S. M. Wise, Analysis of a Darcy-Cahn-Hilliard diffuse interface model for the Hele-Shaw flow and its fully discrete finite element approximation, SIAM Journal on Numerical Analysis, 50 (2012), 1320–1343. [19] J. Gopalakrishnan and G. Kanschat, A multilevel discontinuous Galerkin method, Numerische Mathematik, 95 (2003), 527–550. [20] W. Hackbusch, “Multi-Grid Methods and Applications,” Springer Series in Computational Mathematics. Springer, Berlin, 2010. [21] M. R. Hanisch, Multigrid preconditioning for the biharmonic Dirichlet problem, SIAM Journal on Numerical Analysis, 30 (1993), 184–214. [22] Z. Hu, S. M. Wise, C. Wang and J. S. Lowengrub, Stable and efficient finite-difference nonlinear-multigrid schemes for the phase field crystal equation, Journal of Computational Physics, 228 (2009), 5323–5339. [23] O. A. Karakashian and W. N. Jureidini, A nonconforming finite element method for the stationary Navier-Stokes equations, SIAM Journal on Numerical Analysis, 35 (1998), 93– 120. [24] D. Kay, V. Styles and E. S¨ uli, Discontinuous Galerkin finite element approximation of the Cahn-Hilliard equation with convection, SIAM Journal on Numerical Analysis, 47 (2009), 2660–2685. [25] T. Ohta and K. Kawasaki, Equilibrium morphology of block copolymer melts, Macromolecules, 19 (1986), 2621–2632. [26] P. Percell and M. F. Wheeler, A local residual finite element procedure for elliptic equations, SIAM Journal on Numerical Analysis, 15 (1978), 705–714. [27] U. Trottenberg, C. W. Oosterlee and A. Sch¨ uller, “Multigrid,” Academic Press, New York, 2005. [28] C. Wang, X. Wang and S. M. Wise, Unconditionally stable schemes for equations of thin film epitaxy, Discrete and Continuous Dynamical Systems - Series A (DCDS-A), 28 (2010), 405–423. [29] C. Wang and S. M. Wise, An energy stable and convergent finite-difference scheme for the modified phase field crystal equation, SIAM Journal on Numerical Analysis, 49 (2011), 945– 969. [30] G. N. Wells, E. Kuhl and K. Garikipati, A discontinuous Galerkin method for the CahnHilliard equation, Journal of Computational Physics, 218 (2006), 860–877. [31] M. F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM Journal on Numerical Analysis, 15 (1978), 152–161. [32] S. M. Wise, Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn-Hilliard-Hele-Shaw system of equations, Journal of Scientific Computing, 44 (2010), 38–68. [33] S. M. Wise, C. Wang and J. S. Lowengrub, An energy-stable and convergent finite-difference scheme for the phase field crystal equation, SIAM Journal on Numerical Analysis, 47 (2009), 2269–2288. [34] Y. Xia, Y. Xu and C. W. Shu, Local discontinuous Galerkin methods for the Cahn-Hilliard type equations, Journal of Computational Physics, 227 (2007), 472–491.

Received June 2013; revised August 2013. E-mail address: [email protected] E-mail address: [email protected] E-mail address: [email protected]