A Hybrid High-Order method for the Cahn–Hilliard problem in mixed formThis work was partially supported by Saint-Gobain Recherche (contract UM 150095). D. Di Pietro also acknowledges the partial support of Agence Nationale de la Recherche project HHOMM (ANR-15-CE40-0005).

A Hybrid High-Order method for the Cahn–Hilliard problem in mixed form1

Abstract

In this work we develop a fully implicit Hybrid High-Order algorithm for the Cahn–Hilliard problem in mixed form. The space discretization hinges on local reconstruction operators from hybrid polynomial unknowns at elements and faces. The proposed method has several advantageous features: (i) It supports fairly general meshes possibly containing polyhedral elements and nonmatching interfaces; (ii) it allows arbitrary approximation orders; (iii) it has a moderate computational cost thanks to the possibility of locally eliminating element-based unknowns by static condensation. We perform a detailed stability and convergence study, proving optimal convergence rates in energy-like norms. Numerical validation is also provided using some of the most common tests in the literature.
2010 Mathematics Subject Classification: 65N08, 65N30, 65N12
Keywords: Hybrid High-Order, Cahn–Hilliard equation, phase separation, mixed formulation, discrete functional analysis, polyhedral meshes

\headers

A HHO method for the Cahn–Hilliard problemF. Chave, D. A. Di Pietro, F. Marche, F. Pigeonneau

1 Introduction

Let , , denote a bounded connected convex polyhedral domain with boundary and outward normal , and let . The Cahn–Hilliard problem, originally introduced in [11, 10] to model phase separation in a binary alloy, consists in finding the order-parameter and chemical potential such that

(1a)
(1b)
(1c)
(1d)

where such that on denotes the initial datum, the interface parameter (usually taking small values), and the free-energy such that

(2)

Relevant extensions of problem (1) (not considered here) include the introduction of a flow which requires, in particular, to add a convective term in (1a); cf., e.g., [29, 5, 7, 8, 31, 30].

The discretization of the Cahn–Hilliard equation (1) has been considered in several works. Different aspects of standard finite element schemes have been studied, e.g., in [22, 21, 14]; cf. also the references therein. Mixed finite elements are considered in [24]. In [35], the authors study a nonconforming method based on shape functions for the fourth-order primal problem obtained by plugging (1b) into (1a). Discontinuous Galerkin (dG) methods have also received extensive attention. We can cite here [36], where a local dG method is proposed for a Cahn–Hilliard system modelling multi-component alloys, and a stability analysis is carried out; [23], where optimal error estimates are proved for a dG discretization of the Cahn–Hilliard problem in primal form; [30], which contains optimal error estimates for a dG method based on the mixed formulation of the problem including a convection term;  [26], where a multi-grid approach is proposed for the solution of the systems of algebraic equations arising from a dG discretization of the Cahn–Hilliard equation. In all of the above references, standard meshes are considered. General polygonal meshes in dimension , on the other hand, are supported by the recently proposed -conforming Virtual Element (VE) method of [4] for the problem in primal formulation; cf. also [6] for VE methods with arbitrary regularity. Therein, the convergence analysis is carried out under the assumption that the discrete order-parameter satisfies a -like a priori bound.

In this work, we develop and analyze a fully implicit Hybrid High-Order (HHO) algorithm for problem (1) where the space discretization is based on the HHO( variation proposed in [12] of the method of [19]. The method hinges on hybrid degrees of freedom (DOFs) located at mesh elements and faces that are polynomials of degree and , respectively. The nonlinear term in (1b) is discretized by means of element unknowns only. For the second-order diffusive operators in (1a) and (1b), on the other hand, we rely on two key ingredients devised locally inside each element: (i) A potential reconstruction obtained from the solution of (small) Neumann problems and (ii) a stabilization term penalizing the lowest-order part of the difference between element- and face-based unknowns. See also [13, 34, 33] for related methods for second-order linear diffusion operators, each displaying a set of distinctive features. The global discrete problem is then obtained by a standard element-by-element assembly procedure. When using a first-order (Newton-like) algorithm to solve the resulting system of nonlinear algebraic equations, element-based unknowns can be statically condensed. As a result, the only globally coupled unknowns in the linear subproblems are discontinuous polynomials of degree on the mesh skeleton for both the order-parameter and the chemical potential. With a backward Euler scheme to march in time, the -like error on the order-parameter and the -like error on the chemical potential are proved to optimally converge as (with and denoting, respectively, the spatial and temporal mesh sizes) provided the solution has sufficient regularity.

The proposed method has several advantageous features: (i) It supports general meshes possibly including polyhedral elements and nonmatching interfaces (resulting, e.g., from nonconforming mesh refinement); (ii) it allows one to increase the spatial approximation order to accelerate convergence in the presence of (locally) regular solutions; (iii) it is (relatively) inexpensive. When , e.g., the number of globally coupled spatial unknowns for our method scales as (with denoting the number of mesh faces) as opposed to (with denoting the number of mesh elements) for a mixed dG method delivering the same order of convergence (i.e., based on broken polynomials of degree ). Additionally, thanks to the underlying fully discontinuous polynomial spaces, the proposed method can accomodate abrupt variations of the unknowns in the vicinity of the interface between phases.

Our analysis adapts the techniques originally developed in [30] in the context of dG methods. Therein, the treatment of the nonlinear term in (1b) hinges on -in-time a priori estimates for various norms and seminorms of the discrete order-parameter. Instrumental in proving these estimates are discrete functional analysis results, including discrete versions of Agmon’s and Gagliardo–Nirenberg–Poincaré’s inequalities for broken polynomial functions on quasi-uniform matching simplicial meshes. Adapting these tools to hybrid polynomial spaces on general meshes entails several new ideas. A first key point consists in defining appropriate discrete counterparts of the Laplace and Green’s operators. To this end, we rely on a suitably tailored -like hybrid inner product which guarantees stability estimates for the former and optimal approximation properties for the latter. Another key point consists in replacing the standard nodal interpolator used in the proofs of [30, Lemmas 2.2 and 2.3] by the -orthogonal projector which, unlike the former, is naturally defined for meshes containing polyhedral elements. We show that this replacement is possible thanks to the -stability and approximation properties of the -orthogonal projector on broken polynomial spaces recently presented in a unified setting in [15]; cf. also the references therein for previous results on this subject.

The material is organized as follows: In Section 2 we introduce the notation for space and time meshes and recall some key results on broken polynomial spaces; in Section 3 we introduce hybrid polynomial spaces and local reconstructions, and state the discrete problem; in Section 4 we carry out the stability analysis of the method, while the convergence analysis is detailed in Section 5; Section 6 contains an extensive numerical validation of the proposed algorithm; finally, in Appendix A we give proofs of the discrete functional analysis results used to derive stability bounds and error estimates.

2 Discrete setting

In this section we introduce the discrete setting and recall some basic results on broken polynomial spaces.

2.1 Space and time meshes

We recall here the notion of admissible spatial mesh sequence from [17, Chapter 1]. For the sake of simplicity, we will systematically use the term polyhedral also when . Denote by a countable set of spatial meshsizes having as its unique accumulation point. We consider -refined mesh sequences where, for all , is a finite collection of nonempty disjoint open polyhedral elements of boundary such that and with standing for the diameter of the element .

A face is defined as a planar closed connected subset of with positive -dimensional Hausdorff measure and such that (i) either there exist such that and is called an interface or (ii) there exists such that and is called a boundary face. Mesh faces are collected in the set , and the diameter of a face is denoted by . For all , denotes the set of faces lying on and, for all , is the unit normal to pointing out of . Symmetrically, for all , we denote by the set of one (if ) or two (if ) elements sharing .

{assumption}

[Admissible spatial mesh sequence] We assume that, for all , admits a matching simplicial submesh and there exists a real number independent of such that, for all , the following properties hold: (i) Shape regularity:For all simplex of diameter and inradius , ; (ii) contact-regularity:For all , and all such that , .

To discretize in time, we consider a uniform partition of the time interval with , and for all (the analysis can be adapted to nonuniform partitions). For any sufficiently regular function of time taking values in a vector space , we denote by its value at discrete time , and we introduce the backward differencing operator such that, for all ,

(3)

In what follows, we often abbreviate by the inequality with and positive real numbers and generic constant independent of both the meshsize and the time step (named constants are used in the statements for the sake of easy consultation). Also, for a subset , we denote by and the usual -inner product and norm, with the convention that we omit the index if . The same notation is used for the vector-valued space .

2.2 Basic results on broken polynomial spaces

The proposed method is based on local polynomial spaces on mesh elements and faces. Let an integer be fixed. Let be a subset of , the affine space spanned by , its dimension, and assume that has a non-empty interior in . We denote by the space spanned by -variate polynomials on of total degree , and by the -orthogonal projector onto this space. In the following sections, the set will represent a mesh element or face. The space of broken polynomial functions on of degree is denoted by , and is the corresponding -orthogonal projector.

We next recall some functional analysis results on polynomial spaces. The following discrete trace and inverse inequalities are proved in [17, Chapter 1] (cf. in particular Lemmas 1.44 and 1.46): There is independent of such that, for all , and all ,

(4)

and

(5)

We will also need the following local direct and reverse Lebesgue embeddings (cf. [15, Lemma 5.1]): There is independent of such that, for all , all ,

(6)

The proof of the following results for the local -orthogonal projector can be found in [15, Appendix A.2]. For an open set of , and , we define the seminorm as follows: For all ,

where and . For , we recover the usual Lebesgue spaces . The -orthogonal projector is -stable and has optimal -approximation properties: There is independent of such that, for all , all , all , and all , it holds,

(7)

and, for all ,

(8)

where denotes the set of functions that belong to for all . Finally, there is independent of such that it holds, for all ,

(9)

In the proofs of Lemmas 1 and 4 below, we will make use of the following global inverse inequalities, which require mesh quasi-uniformity.

Proposition 1 (Global inverse inequalities for Lebesgue norms of broken polynomials)

In addition to Assumption 2.1, we assume that the mesh is quasi-uniform, i.e.,

(10)

Then, for all polynomial degree and all , it holds

(11)

with real number independent of .

Proof

Let . We start by proving that, for all ,

(12)

which corresponds to (11) with . By the local reverse Lebesgue embeddings (6), there is independent of such that

where we have used the mesh quasi-uniformity assumption (10) to conclude. Inequality (12) follows observing that . Let us now turn to the case . We have

where the conclusion follows using (12).

3 The Hybrid High-Order method

In this section we define hybrid spaces and state the discrete problem.

3.1 Hybrid spaces

The discretization of the diffusion operator hinges on the HHO method of [12] using polynomials of degree inside elements and on mesh faces (cf. Remark 6 for further insight on this choice). The global discrete space is defined as

(13)

The restriction of to an element is denoted by . For a generic collection of DOFs in , we use the underlined notation and, for all , we denote by its restriction to . Also, to keep the notation compact, we denote by (no underline) the function in such that

In what follows, we will also need the zero-average subspace

The interpolator is such that, for all ,

(14)

We define on the seminorm such that

(15)

where denotes the usual broken gradient on and the stabilization bilinear form on is such that

(16)

Using the stability and approximation properties of the -orthogonal projector expressed by (7)–(8), one can prove that is -stable:

(17)

The following Friedrichs’ inequalities can be proved using the arguments of [15, Lemma 7.2], where element DOFs of degree are considered (cf. also [9, 16] for related results using dG norms): For all if , all if ,

(18)

The case corresponds to Poincaré’s inequality. Finally, to close this section, we prove that defines a norm on .

Proposition 2 (Norm )

The map defines a norm on .

Proof

We only have to show that . By (18), , i.e., for all . Plugging this result into the definition (15) of , we get , which implies in turn for all .

3.2 Diffusive bilinear form and discrete problem

For all , we define the potential reconstruction operator such that, for all , is the unique solution of the following Neumann problem:

(19)

with closure condition . It can be proved that, for all , denoting by the restriction of the reduction map defined by (14) to ,

(20)

which expresses the fact that is the elliptic projector onto (and, as such, has optimal approximation properties in ). The diffusive bilinear form on is obtained by element-wise assembly setting

(21)

with stabilization bilinear form defined by (16). Denoting by the seminorm defined by on , a straightforward adaptation of the arguments used in [19, Lemma 4] shows that

(22)

which expresses the coercivity and boundedness of . Additionally, following the arguments in [19, Theorem 8], one can easily prove that the bilinear form enjoys the following consistency property: For all , , such that on ,

(23)
Remark 1 (Consistency of )

For sufficiently regular solutions (i.e., when ), equation (23) shows that the consistency error scales as . This is a consequence of the fact that both the potential reconstruction (cf. (19)) and the stabilization bilinear form (cf. (16)) are consistent for exact solutions that are polynomials of degree inside each element. In particular, a key point in is to penalize instead of . A similar stabilization bilinear form had been independently suggested in the context of Hybridizable Discontinuous Galerkin methods in [32, Remark 1.2.4].

The discrete problem reads: For all , find such that

(24a)
(24b)

and solves

(25)

We note, in passing, that the face DOFs in are not needed to initialize the algorithm.

Remark 2 (Static condensation)

Problem (24) is a system of nonlinear algebraic equations, which can be solved using an iterative algorithm. When first order (Newton-like) algorithms are used, element-based DOFs can be locally eliminated at each iteration by a standard static condensation procedure.

4 Stability analysis

In this section we establish some uniform a priori bounds on the discrete solution. To this end, we need a discrete counterpart of Agmon’s inequality; cf. [3, Lemma 13.2] and also [1, Theorem 3]. We define on the following -like inner product:

(26)

and denote by and the norm and seminorm corresponding to the bilinear forms and , respectively. For further insight on the role of , cf. Remark 7. We introduce the discrete Laplace operator such that, for all ,

(27)

and we denote by (no underline) the broken polynomial function in obtained from element DOFs in .

Remark 3 (Restriction of to )

Whenever , . To prove it, it suffices to take in (27) (with characteristic function of ), and observe that the left-hand side satisfies while, by definition (21) of the bilinear form , the right-hand side vanishes. In what follows, we keep the same notation for the (bijective) restriction of to .

The following result, valid for , will be proved in Appendix A.

Lemma 1 (Discrete Agmon’s inequality)

Assume mesh quasi-uniformity (10). Then, it holds with real number independent of ,

(28)

We also recall the following discrete Gronwall’s inequality (cf. [28, Lemma 5.1]).

Lemma 2 (Discrete Gronwall’s inequality)

Let two reals be given, and, for integers , let , , and denote nonnegative real numbers such that

Then, if for all , letting , it holds

(29)

We are now ready to prove the a priori bounds.

Lemma 3 (Uniform a priori bounds)

Under the assumptions of Lemma 1, and further assuming that for a given real number independent of and of (but depending on ) and sufficiently small, there is a real number independent of and such that

Proof

The proof is split into several steps.

  1. We start by proving that

    (30)

    Subtracting (24b) with from (24a) with , and using the fact that, for all , , it is inferred, for all , that

    (31)

    Notice that for all by definition (2) of . Making in (24a) and using the Cauchy–Schwarz and Young’s inequalities, we infer that

    (32)

    Additionally, recalling the following formula for the backward Euler scheme:

    (33)

    it holds

    (34)

    Plugging (32) and (34) into (31), we obtain

    Provided , the bound (30) follows summing the above inequality over , and using the fact that . To prove this bound, observe that

    where we have used the definition (2) of the free-energy in the first line followed by the discrete Friedrichs’ inequality with in the second line and the first bound on the initial datum in (46) below to conclude.

  2. We next prove that

    (35)

    The discrete Agmon’s inequality (28) followed by the first inequality in (22) yields