Local Discontinuous Galerkin Methods for the Functionalized Cahn–Hilliard Equation

Abstract

In this paper, we develop a local discontinuous Galerkin (LDG) method for the sixth order nonlinear functionalized Cahn–Hilliard (FCH) equation. We address the accuracy and stability issues from simulating high order stiff equations in phase-field modeling. Within the LDG framework, various boundary conditions associated with the background physics can be naturally implemented. We prove the energy stability of the LDG method for the general nonlinear case. A semi-implicit time marching method is applied to remove the severe time step restriction ($\Delta t \sim O(\Delta x^6)$) for explicit methods. The $h-p$ adaptive capability of the LDG method allows for capturing the interfacial layers and the complicated geometric structures of the solution with high resolution. To enhance the efficiency of the proposed approach, the multigrid (MG) method is used to solve the system of linear equations resulting from the semi-implicit temporal integration at each time step. We show numerically that the MG solver has mesh-independent convergence rates. Numerical simulation results for the FCH equation in two and three dimensions are provided to illustrate that the combination of the LDG method for spatial approximation, semi-implicit temporal integration with the MG solver provides a practical and efficient approach when solving this family of problems.

Appendices

Appendix 1: The Linear MG Solver

We start by introducing the basic notations to describe the general setting of a two-level or bi-grid method. Together with the family of partitions $\{\mathcal {T}_h\}_{h>0}$ used for the LDG discretization, we consider a coarse family of mesh partitions, $\{\mathcal {T}_H\}_{H>0}$ with $H>h$ and satisfying the basic assumption $\mathcal {T}_H \subset \mathcal {T}_h$. One can think $H=2h$, since in many circumstances it will already be coarse enough. Associated to the coarse mesh partition we have the corresponding finite element space $V_H$ which is defined as

$$\begin{aligned} V_H=\{v \in \, L^{2}(D): v|_D\in \mathcal {P}^k(D);\quad \forall \, D\in \mathcal {T}_H\}. \end{aligned}$$

(8.1)

Throughout the whole description we assume the polynomial degree $k$ is fixed.

To link functions in both spaces, we define the prolongation and restriction operators. The prolongation operator $P_{hH}:V_H \vee V_h$ is defined as the natural inclusion. The restriction operator $R_{Hh} : V_h \vee V_H$ is defined as the transpose of $P_{hH}$ with respect to the standard $L^{2}$-inner product. That is, it is obtained by solving:

$$\begin{aligned} \sum _{D\in \mathcal {T}_H}\int _{D} R_{Hh}(u_h) v_{H} dx=\sum _{D\in \mathcal {T}_h} \int _{D} u_h P_{hH}(v_{H}) dx,\quad \forall v_{H} \in V_H. \end{aligned}$$

(8.2)

We now denote $S_h$ as a general relaxation or smoothing operator. We later specify and study several choices. The basic property that $S_h$ should have is to damp the high frequencies of the approximate solution and smooth the error. The coarse solver is defined by $A_{H}= R_{Hh}A_hP_{hH}$. The coarse grid correction step is to reduce the smooth components of the error that can not be reduced by the smoother.

The linear scheme requires the solution at each time step of algebraic non-symmetric linear systems i.e.

$$\begin{aligned} A_hu_h=f_h. \end{aligned}$$

(8.3)

The MG method is used to solve the system and the main points of the algorithm is the bi-grid cycle. Following [2], we can formulate the bi-grid cycle as follows:

Algorithm 1: Two-grid Cycle

Starting with an initial approximation, say $u_{PRE}^0$:

• 1. Pre-relaxation apply $\nu _1$ pre-relaxation sweeps: for $m=1\ldots , \nu _1$, solve

  $$\begin{aligned} u_{PRE}^{m}=u_{PRE}^{m-1}+S_h(f_h-A_h u_{PRE}^{m-1}), \end{aligned}$$

• 2. Coarse-grid correction update the solution $u_{PRE}^{\nu _1}$ by a coarse-grid correction step, solve the problem once on coarse grid

  $$\begin{aligned} A_{H}v_{H}=R_{Hh}(f_h-A_h u_{PRE}^{\nu _1}), \end{aligned}$$

  and set $u_{CG}=u_{PRE}^{\nu _1}+ P_{hH}(v_{H})$.

• 3. Post-relaxation starting with $u_{CG}$, apply $\nu _2$ post-relaxation sweeps, that is, set $u_{POST}^{0}=u_{CG}$ and for $m=1\ldots , \nu _2$, solve

  $$\begin{aligned} u_{POST}^{m}=u_{POST}^{m-1}+S_h(f_h-A_hu_{POST}^{m-1}). \end{aligned}$$

The integers $\nu _1$ and $\nu _2$ are parameters in the scheme that control the number of relaxation sweeps before and after visiting the coarse grid. $\nu _1$ and $\nu _2$ are called the number of pre- and post- relaxations, respectively.

Solving the coarse grid problem at the second step of the above algorithm could be done again with the two-level algorithm. Hence, the V-cycle multi-level algorithm in terms of the two-level algorithm is defined by applying the two-level algorithm recursively.

Due to the block structure of $A_h$, we focus only on very simple block-relaxation. In particular, we decompose $A_h$ into a strict block-lower, a block-diagonal, and a strict block-upper matrix, i.e.

$$\begin{aligned} A_h=L_h+D_h+U_h. \end{aligned}$$

(8.4)

Table 4 shows some possible choices of the smoother operator.

Table 4 The possible choices of the smoother operator

Appendix 2: The Local Mode Analysis of the Two-Grid Algorithm

In [1], the author considered a general framework for performing the local mode analysis for analyzing the convergence of two-level or bi-grid algorithms, and also provided some quantitative information about the performance and design of the solvers. Although the approach is applied to constant coefficients, linear problems and uniform grids, some general results are established, based on the fact that some local mode analysis can be performed at the matrices. A different treatment has to be given to the part of the matrix associated with interior unknowns and that associated to boundary degrees of freedom. Ignoring the treatment of the boundaries, we now revise some of the results given in [1]. There, the author defined the convergence factor of the two-grid method by

$$\begin{aligned} \lambda =\sup \frac{\Vert u_{POST}^{\nu _2}\Vert }{\Vert u_{PRE}^{0}\Vert } \end{aligned}$$

(9.1)

in some appropriate chosen norm that might depend on the problem. It is also shown, that under the assumptions described above (constant coefficients, linear problems, uniform grids and neglecting the boundary conditions), the convergence factor might be computed in terms of the symbol of the error propagation operator.

For the linear iteration, the error propagation operator is defined as:

$$\begin{aligned} E_h:= I- S_hA_h, \end{aligned}$$

(9.2)

where $I$ is the identity operator in $V_h^{k}$. The spectral radius, or some norm of this operator allows to quantify how the error is reduced at each iteration. If it is less than $1$ we will get a convergent iteration. The smaller it is, the faster is the iteration.

In the case of the two level or two-grid cycle defined in Algorithm 1, the error propagation is:

$$\begin{aligned} E_h^{2grid}:= E_h^{\nu _2} [I- P_{hH}A_{H}^{-1}R_{Hh}A_h]E_h^{\nu _1}. \end{aligned}$$

(9.3)

Following [1], if one could choose the $L^{2}$-norm, denotes by $\widehat{E}_h(\theta )$ the symbol (in the frequency space) of the error propagation operator, from Parseval’s identity it is formally obtained

$$\begin{aligned} \lambda =\sup _{\theta \ne 0} \Vert \widehat{E}_h(\theta )\Vert . \end{aligned}$$

(9.4)

While for symmetric problems, the estimation of the spectral radius of $E_h$ could be reduced to the computation of its largest eigenvalue, in the present situation, since $A_h$ is non-symmetric and also $S_h$, one can not guarantee that their spectral information contain the relevant information.

To compute the spectral radius, a possible way is to compute the first singular value of $E_h$. In particular, one can define the asymptotic convergence factor (see [1]) as:

$$\begin{aligned} \lambda _{asymp}=\sup _{\theta \ne 0} \sigma _1(\widehat{E}_h(\theta )), \end{aligned}$$

(9.5)

where $\sigma _1$ is the spectral radius of $E$, (i.e. largest absolute eigenvalue). On the other hand, a more restrictive way of ensuring that $\Vert E_h^{2grid}\Vert <1$ would be to study the norms of each of the terms in the product. From the definition, we see that an essential requirement is to guarantee that the smoother has norm strictly less than $1$.