Abstract
In this paper, we develop and analyze a fast solver for the system of algebraic equations arising from the local discontinuous Galerkin (LDG) discretization and implicit time marching methods to the Cahn–Hilliard (CH) equations with constant and degenerate mobility. Explicit time marching methods for the CH equation will require severe time step restriction \((\varDelta t \sim O(\varDelta x^4))\), so implicit methods are used to remove time step restriction. Implicit methods will result in large system of algebraic equations and a fast solver is essential. The multigrid (MG) method is used to solve the algebraic equations efficiently. The Local Mode Analysis method is used to analyze the convergence behavior of the linear MG method. The discrete energy stability for the CH equations with a special homogeneous free energy density \(\Psi (u)=\frac{1}{4}(1-u^2)^2\) is proved based on the convex splitting method. We show that the number of iterations is independent of the problem size. Numerical results for one-dimensional, two-dimensional and three-dimensional cases are given to illustrate the efficiency of the methods. We numerically show the optimal complexity of the MG solver for \(\mathcal{P }^1\) element. For \(\mathcal{P }^2\) approximation, the optimal or sub-optimal complexity of the MG solver are numerically shown.

















Similar content being viewed by others
References
Alexander, R.: Diagonally implicit Runge–Kutta methods for stiff O.D.E’.S. SIAM J. Numer. Anal. 14, 1006–1021 (1977)
Barrett, J.W., Blowey, J.F.: Finite element approximation of a model for phase separation of a multi-component alloy with non-smooth free energy. Numer. Math. 77, 1–34 (1997)
Barrett, J.W., Blowey, J.F., Garcke, H.: Finite element approximation of the Cahn–Hilliard equation with degenerate mobility. SIAM J. Numer. Anal. 37, 286–318 (1999)
Barrett, J.W., Blowey, J.F., Garcke, H.: On fully practical finite element approximations of degenerate Cahn–Hilliard systems. M2ANMath. Model. Numer. Anal. 35, 713–748 (2001)
Bassi, F., Ghidoni, A., Rebay, S., Tesini, P.: High-order accurate \(p\)-multigrid discontinuous Galerkin solution of the Euler equations. Int. J. Numer. Methods Fluids 60, 847–865 (2008)
Bassi, F., Ghidoni, A., Rebay, S.: Optimal Runge–Kutta smoothers for the \(p\)-multigrid discontinuous Galerkin solution of the 1D Euler equations. J. Comput. Phys. 230, 4153–4175 (2011)
Blowey, J.F., Elliott, C.M.: The Cahn–Hilliard gradient theory for phase separation with nonsmooth free energy. II. Numerical analysis. Eur. J. Appl. Math. 3, 147–179 (1992)
Blowey, J.F., Copetti, M.I.M., Elliott, C.M.: Numerical analysis of a model for phase separation of multi-component alloy. IMA J. Numer. Anal. 16, 111–139 (1996)
Brandt, A.: Multigrid techniques: 1984 guide with applications to fluid dynamics. GMD-Studien [GMD Studies], 85. Gesellschaft für Mathematik und Datenverarbeitung mbH, St. Augustin (1984)
Brandt, A.: Rigorous quantitative analysis of multigrid. I. Constant coefficients two-level cycle with \(L_2\)-norm. SIAM J. Numer. Anal. 31, 1695–730 (1994)
Choo, S.M., Lee, Y.J.: A discontinuous Galerkin method for the Cahn–Hilliard equation. J. Appl. Math. Comput. 18, 113–126 (2005)
Cockburn, B., Shu, C.-W.: The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal. 35, 2440–2463 (1998)
Copetti, M.I.M., Elliott, C.M.: Numerical analysis of the Cahn–Hilliard equation with logarithmic free energy. Numer. Math. 63, 39–65 (1992)
Elliott, C.M., French, D.A.: A nonconforming finite-element method for the two-dimensional Cahn–Hilliard equation. SIAM J. Numer. Anal. 26, 884–903 (1989)
Elliott, C.M., French, D.A., Milner, F.A.: A second order splitting method for the Cahn–Hilliard equation. Numer. Math. 54, 575–590 (1989)
Eyre, D.J.: Systems of Cahn–Hilliard equations. SIAM J. Appl. Math. 53, 1686–1712 (1993)
Feng, X., Prohl, A.: Analysis of a fully discrete finite element method for the phase field model and approximation of its sharp interface limits. Math. Comput. 73, 541–567 (2004)
Feng, X.B., Karakashian, O.A.: Fully discrete dynamic mesh discontinuous Galerkin methods for the Cahn–Hilliard equation of phase transition. Math. Comput. 76, 1093–1117 (2007)
Furihata, D.: A stable and conservative finite difference scheme for the Cahn–Hilliard equation. Numer. Math. 87, 675–699 (2001)
Hector, Gomez B., Calo, Victor M., Hughes Thomas, J.R.: Isogeometric analysis of the Cahn–Hilliard phase-field model. Comput. Methods Appl. Mech. Eng. 197, 4333–4352 (2008)
Guo, R., Xu, Y., Blanca Ayuso de Dios.: Multigrid methods of discontinuous Galerkin discretization for linear time-dependent fourth-order equations, Preprint
Kay, D., Welford, R.: A multigrid finite element solver for the Cahn–Hilliard equation. J. Comput. Phys. 212, 288–304 (2006)
Kim, J., Kang, K., Lowengrub, J.: Conservative multigrid methods for Cahn–Hilliard fluids. J. Comput. Phys. 193, 511–543 (2004)
Kim, J., Kang, K., Lowengrub, J.: Conservative multigrid methods for ternary Cahn–Hilliard systems. Commun. Math. Sci. 2, 53–77 (2004)
Klaij, C.M., van Raalte, M.H., van der Ven, H., van der Vegt, J.J.W.: \(h\)-Multigrid for space-time discontinuous Galerkin discretizations of the compressible Navier–Stokes equations. J. Comput. Phys. 227, 1024–1045 (2007)
Reed, W., Hill, T.: Triangular mesh methods for the neutrontransport equation, La-ur-73-479, Los Alamos Scientific Laboratory (1973)
Shahbazi, K., Mavriplis, D.J., Burgess, N.K.: Multigrid algorithms for high-order discontinuous Galerkin discretizations of the compressible Navier–Stokes equations. J. Comput. Phys. 228, 7917–7940 (2009)
Sun, Z.Z.: A second-order accurate linearized difference scheme for the twodimensional Cahn–Hilliard equation. Math. Comput. 64, 1463–1471 (1995)
Trottenberg, U., Oosterlee, C., Schller, A.: Multigrid. Academic Press, New York (2005)
van der Vegt, J.J.W., Rhebergen, S.: \(hp\)-multigrid as smoother algorithm for higher order discontinuous Galerkin discretizations of advection dominated flows. Part I. Multilevel analysis. J. Comput. Phys. 231, 7537–7563 (2012)
van der Vegt, J.J.W., Rhebergen, S.: \(hp\)-multigrid as smoother algorithm for higher order discontinuous Galerkin discretizations of advection dominated flows. Part II. Optimization of the Runge–Kutta smoother. J. Comput. Phys. 231, 7564–7583 (2012)
Wells, G.N., Kuhl, E., Garikipati, K.: A discontinuous Galerkin method for the Cahn–Hilliard equation. J. Comput. Phys. 218, 860–877 (2006)
Wise, S.M.: Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn–Hilliard–Hele–Shaw system of equations. J. Sci. Comput 44, 38–68 (2010)
Xia, Y., Xu, Y., Shu, C.-W.: Local discontinuous Galerkin methods for the Cahn–Hilliard type equations. J. Comput. Phys. 227, 472–491 (2007)
Xia, Y., Xu, Y., Shu, C.-W.: Efficient time discretization for local discontinuous Galerkin methods. Discret. Contin. Dyn. Syst. Ser. B 8, 677–693 (2007)
Xu, Y., Shu, C.-W.: Local discontinuous Galerkin methods for high-order time-dependent partial differential equations. Commun. Comput. Phys. 7, 1–46 (2010)
Author information
Authors and Affiliations
Corresponding author
Additional information
Research supported by NSFC Grant No. 10971211, No. 11031007, FANEDD No. 200916, NCET No. 09-0922, Fok Ying Tung Education Foundation No. 131003.
Appendices
Appendix 1: The Linear MG Solver
The linear scheme requires the solution at each time step of algebraic non-symmetric linear systems i.e.
The MG method is used to solve the system and the main points of the algorithm is the bi-grid cycle. Following [9], 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-gird 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.
Table 12 shows some possible choices of the smoother operator.
Appendix 2: The Nonlinear FAS MG Solver
The nonlinear scheme (3.4) requires the solution at each time step of three algebraic nonlinear systems, i.e.
and the nonlinear FAS MG method is introduced to solve the system. Following [29], we illustrate the FAS bi-grid cycle as follows:
Algorithm 2: Two-grid Cycle Starting with an initial approximation, say \(u_h^m \):
-
1.
Pre-relaxation: apply \(\nu _1\) pre-relaxation sweeps to (8.1) with \(u_h^m \) the initial approximation and obtain \(\bar{u}_h^m \).
-
2.
Coarse-gird correction: compute the coarse level initial iterate
$$\begin{aligned} \bar{u}_H^m=R_{Hh}\bar{u}_h^m, \end{aligned}$$(8.2)update the solution \(\bar{u}_h^m \) by a coarse-grid correction step, solve the problem once on coarse grid
$$\begin{aligned} N_H(v_H)=N_H(\bar{u}_H^m)+R_{Hh}(f_h-N_h(\bar{u}_h^m)) \end{aligned}$$(8.3)and set \(u_h^{CG}=\bar{u}_h^m+P_{hH}(v_H-\bar{u}_H^m)\).
-
3.
Post-relaxation: starting with \(u_h^{CG}\), apply \(\nu _2\) post-relaxation sweeps to (8.1) and obtain \(u_h^{m+1} \).
Solving the coarse grid problem at the second step of the above algorithm could be done again with the two-level algorithm. Hence, the nonlinear FAS MG algorithm in terms of the two-level algorithm is defined by applying the two-level algorithm recursively.
Appendix 3: The Local Mode Analysis
In [10], 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 [10]. There, the author defined the convergence factor of the two-grid method by
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 condition), 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:
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:
Following [10], 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
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 [10]) as:
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\).
Rights and permissions
About this article
Cite this article
Guo, R., Xu, Y. Efficient Solvers of Discontinuous Galerkin Discretization for the Cahn–Hilliard Equations. J Sci Comput 58, 380–408 (2014). https://doi.org/10.1007/s10915-013-9738-4
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-013-9738-4
Keywords
- Cahn–Hilliard equation
- Local discontinuous Galerkin method
- Convex splitting method
- Multigrid algorithm
- FAS multigrid
- Additive Runge–Kutta
- Diagonally implicit Runge–Kutta
- Local mode analysis