A Discontinuous Galerkin Scheme for the Cahn-Hilliard Equations with Discrete Maximum Principle for Arbitrary Polynomial Order

Abstract: We propose a structure-preserving discontinuous Galerkin scheme for the Cahn--Hilliard equations with degenerate mobility based on the Symmetric Weighted Interior Penalty formulation. By evaluating the mobility at cell averages rather than as a piecewise polynomial, the proposed scheme preserves strict degeneracy and yields a coercivity constant that is independent of the mobility, removing the need for regularisation. Moreover, we establish existence of discrete solutions even with degeneracy via a Leray--Schauder fixed-point argument, and show that the scheme satisfies a provable discrete maximum principle at arbitrary polynomial order $p$ when combined with the Zhang--Shu scaling limiter for $p > 0$ and from the scheme alone for $p = 0$. Mass conservation and energy dissipation are established for the unlimited scheme; for the limited variant, we discuss observed energy dissipation for $p \geq 1$ and potential theoretical solutions. Numerical experiments confirm optimal convergence rates of order $p+1$ in $L^2$ and validate structure-preserving properties with numerical results.

• The paper presents the SWIPDP-L DG scheme that unconditionally preserves a discrete maximum principle for Cahn–Hilliard equations with degenerate mobility.
• It employs cell-average evaluation and the Zhang–Shu scaling limiter to enforce strict bound preservation while conserving mass.
• Numerical experiments verify optimal convergence rates, robust energy dissipation, and exact mass conservation across various polynomial orders.

Discontinuous Galerkin Methods with Discrete Maximum Principle for Cahn–Hilliard Equations

Introduction

The paper introduces a discontinuous Galerkin (DG) discretization for the Cahn–Hilliard (CH) equations with degenerate mobility. The central design aims are unconditional preservation of a discrete maximum principle (DMP) at arbitrary polynomial order, structure preservation (including mass conservation and energy dissipation at the discrete level), and robust handling of degenerate nonlinearities without recourse to regularization. The proposed approach, termed the SWIP Diffusive Projection with Limiter (SWIPDP-L), employs the Symmetric Weighted Interior Penalty (SWIP) method with the Zhang–Shu scaling limiter to enforce strict bounds for the phase-field variable. The methodology enables rigorous proofs of a DMP and energy estimates and is supported by comprehensive numerical experiments.

Mathematical Formulation and Scheme Construction

The Cahn–Hilliard system consists of a fourth-order parabolic PDE for a phase-field $\phi$, formulated as

$$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$

on a bounded Lipschitz domain $\Omega\subset\mathbb{R}^d$ ($d=1,2,3$), with $M(\phi)$ the degenerate mobility (e.g., $M(\phi)=\max\{1-\phi^2, 0\}$), $W$ a double-well potential (either quartic or logarithmic). The chief analytic challenge is that the solution is not a priori bounded for quartic $W$, and classical discretizations (FEM, SIPG) cannot ensure a DMP even at the cell-average level, especially with higher-order polynomials.

The main innovation in the DG scheme is the strict evaluation of the mobility function $M$ at cell averages, resulting in a piecewise-constant and strictly degenerate coefficient. This choice decouples the coercivity of the method from the degeneracy, removing the necessity for artificial regularization in the discrete variational problem. The authors construct the SWIPDP trilinear form, establish discrete coercivity (with constants independent of the mobility), and prove mass conservation and a priori energy bounds at the discrete level.

Discrete Maximum Principle and Limiter Application

Ensuring a DMP at arbitrary order requires not only that cell averages remain bounded, but also that the pointwise values of the numerical solution stay within admissible physical bounds $|\phi|\leq 1$. For piecewise-constant approximation ($$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$0) the scheme itself ensures the maximum principle at the cell-average level, which, by construction, extends to the nodal values. However, for $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$1 the Zhang–Shu scaling limiter is employed after each time step. This limiter rescales local deviations from the cell mean to enforce the $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$2 bound, while crucially preserving mass.

The paper rigorously proves that, under the scheme design, the cell-averaged DMP is maintained at each timestep unconditionally in both mesh size and time step. Application of the limiter then ensures the global $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$3 bound, and energy contraction is numerically observed, even though a formal monotonic energy law for the limited ($$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$4) scheme is not proved.

Figure 1: Energy dissipation and bound preservation over time for random spinodal decomposition, demonstrating robust monotone decay of energy and exact mass conservation.

Existence, Coercivity, and Energy Analysis

A significant contribution is the existence result for the fully discrete degenerate DG scheme, proven by a Leray–Schauder fixed-point argument for the regularized problem and passage to the degenerate limit. The coercivity result demonstrates that the penalty parameter can be chosen independently of the cellwise mobility values, eliminating the need for a lower bound on $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$5 and thus avoiding regularization artifacts.

The paper provides explicit energy estimates, guaranteeing the discrete free energy $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$6 is non-increasing for unlimited schemes, and establishes uniform $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$7 control over the solution. For the limited variant, the energy is still shown to remain bounded, albeit monotonic decay is not theoretically guaranteed.

Numerical Results

Comprehensive experiments explore both the accuracy and structure-preserving properties of the SWIPDP-L scheme. Convergence studies with smooth data confirm optimal rates: order $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$8 in $$\partial_t \phi = Pe^{-1} \nabla \cdot ( M(\phi) \nabla \mu ),\quad \mu = W'(\phi) - Cn^2\Delta\phi,$$9 and order $\Omega\subset\mathbb{R}^d$0 in $\Omega\subset\mathbb{R}^d$1. Simulations of prototypical physical settings—spinodal decomposition and bubble merging—demonstrate strict enforcement of phase-field bounds and robust energy dissipation.

Figure 2: Evolution of maximum bound violation for $\Omega\subset\mathbb{R}^d$2 shows errors at the level of machine precision, indicating practical satisfaction of the DMP throughout the simulation.

Figure 3: Snapshots of phase-field evolution during spinodal decomposition illustrate macroscopic phase separation and the preservation of bounds at all times.

Notably, the DMP is satisfied up to solver and floating-point tolerances; observed violations are always at the level of numerical roundoff. The energy dissipation law is verified empirically for both $\Omega\subset\mathbb{R}^d$3 and $\Omega\subset\mathbb{R}^d$4. For the merging bubbles benchmark, mass conservation is exact, bounds remain fixed at the initial extrema, and there are no spurious oscillations.

Figure 4: Monotone energy dissipation, mass invariance, and sharp, frozen bounds in the coalescence of two initial droplets.

Comparison with Alternative Schemes

The proposed SWIPDP-L method is contrasted with recent alternatives:

• The method of Liu et al. admits DMP for arbitrary $\Omega\subset\mathbb{R}^d$5, but requires a two-stage (DG + constrained minimization) procedure, whereas SWIPDP-L applies a single Zhang–Shu limiter.
• Flux upwinding approaches (Acosta et al.) can guarantee DMP for $\Omega\subset\mathbb{R}^d$6, but do not naturally extend to $\Omega\subset\mathbb{R}^d$7 without sacrificing efficiency.
• Previous projection and interior penalty variants (FEM-L, SIPG-L, SWIP-L) require either empirical parameter tuning or do not provide analytic DMP proofs.

The mobility-free penalty parameter and direct cell-average boundedness provide a simpler and more robust path to structure preservation than these alternatives.

Implications and Future Outlook

The rigorous construction and analysis of the SWIPDP-L scheme for the CH equations is a major step towards unconditionally structure-preserving high-order DG discretization of fourth-order nonlinear PDEs with degenerate diffusion. The implications include:

• Reliable and stable simulation of phase-field models for diffuse interface phenomena in multi-component and multi-physics settings.
• Robustness for large polynomial degrees, facilitating $\Omega\subset\mathbb{R}^d$8-adaptivity and higher accuracy for smooth solutions and sharp interface resolution.
• A framework directly applicable to advective and coupled systems, as established by discrete analysis and numerical verifications.

Practical challenges remain. The high nonlinearity and strict degeneracy (absence of regularization) can induce solver ill-conditioning. The performance of the method compared to two-step or energy-optimization limited approaches remains of practical interest. Extending theoretical energy contraction to the $\Omega\subset\mathbb{R}^d$9 limited case and generalizing to additional physical couplings (Cahn–Hilliard–Navier–Stokes, variable or anisotropic mobilities) are open directions.

Conclusion

The SWIPDP-L method provides a rigorous, efficient, and structure-preserving high-order DG scheme for degenerate Cahn–Hilliard equations. It combines sharp monotonicity and discrete mass conservation with optimal accuracy and verified energy behavior. The method presents a solid platform for structure-preserving simulation of multiphase and interfacial phenomena and motivates further investigation into energy-limited higher-order schemes, numerical solvers for degenerate problems, and applications to coupled and adaptive systems.