🔢 mathematics

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

This paper proposes a structure-preserving discontinuous Galerkin scheme for the Cahn-Hilliard equations with degenerate mobility that ensures strict degeneracy without regularization, guarantees existence of discrete solutions, satisfies a discrete maximum principle at arbitrary polynomial orders via a scaling limiter, and maintains mass conservation and energy dissipation, all validated by optimal convergence rates in numerical experiments.

Original authors: Jimmy Kornelije Gunnarsson, Robert Klöfkorn

Published 2026-04-02

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Jimmy Kornelije Gunnarsson, Robert Klöfkorn

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are watching a drop of oil mix with water. At first, they are separate, but over time, they start to swirl, separate, and form intricate patterns like a marble cake. This process is called phase separation, and in the world of physics and engineering, it's described by a complex set of rules known as the Cahn-Hilliard equations.

However, simulating this on a computer is like trying to paint a masterpiece with a shaky hand. The math is tricky because:

The "No-Go" Zones: The mixture can only exist within a specific range (let's say between -1 and 1). If the computer simulation accidentally calculates a value of 1.5 or -2, the whole picture breaks, and the simulation crashes or produces nonsense.
The "Stiff" Problem: In some areas, the mixture stops moving entirely (like oil that has frozen). Standard computer methods struggle with these "dead zones," often requiring the programmer to cheat by adding a tiny bit of fake movement just to keep the math working. This is called "regularization," but it ruins the accuracy.

The New Solution: A Smart, Self-Correcting Painter

The authors of this paper, Jimmy Gunnarsson and Robert Klöforn, have invented a new digital painting technique (a Discontinuous Galerkin scheme) that solves these problems without cheating.

Here is how their method works, using some everyday analogies:

1. The "Cell Average" Rule (The Neighborhood Watch)

Imagine the computer screen is divided into a grid of tiny neighborhoods (cells). In older methods, the computer tried to guess the exact value of the mixture at every single point inside a neighborhood. If one point got too high, the whole neighborhood got confused.

The new method is smarter. It says: "Let's only look at the average mood of the entire neighborhood."

The Trick: They calculate the "mobility" (how easily the mixture moves) based on this average, not on every single point.
The Benefit: If the average is safe, the whole neighborhood is safe. This allows them to handle the "dead zones" (where the mixture stops moving) perfectly without needing to cheat with "regularization." It's like a neighborhood watch that only cares if the average noise level is too high, rather than punishing every single whisper.

2. The "Safety Limiter" (The Bouncer)

Even with the neighborhood rule, sometimes the computer's math might get a little wild and try to push a value slightly outside the safe range (e.g., 1.0000001 instead of 1).

To fix this, they use a Zhang-Shu Scaling Limiter. Think of this as a strict bouncer at a club.

The computer does its calculation first.
Then, the bouncer checks the result. If the mixture value tries to go above 1 or below -1, the bouncer gently pulls it back to the edge of the door.
The Magic: They proved mathematically that if the neighborhood average is safe, the bouncer can always fix the individual points to keep them safe, no matter how complex the math gets. This works for any level of detail (polynomial order), from simple blocks to highly detailed curves.

3. Keeping the Energy Downhill (The Sliding Hill)

In physics, systems naturally lose energy over time (like a ball rolling down a hill until it stops). A good simulation must respect this; it shouldn't make the ball roll up the hill by accident.

The authors proved that their method respects this "energy downhill" rule.
The Catch: For the simplest version of their method, they proved this mathematically. For the more complex, detailed versions, they proved the "bouncer" keeps the values safe, and their computer experiments show the energy still goes downhill, even if a perfect mathematical proof for that specific part is still being worked out.

Why Does This Matter?

No More Cheating: You don't need to add fake "regularization" to make the math work. The simulation is honest about where the mixture stops moving.
Safety First: The simulation will never crash because the values go out of bounds. The "bouncer" ensures the mixture stays within the physical limits (-1 to 1).
High Detail: You can use this method to get very high-resolution, detailed pictures of the mixing process without losing stability.

The Bottom Line

The authors have built a digital simulation engine for mixing fluids that is robust, honest, and safe. It uses a "neighborhood average" strategy to handle difficult physics and a "bouncer" to ensure the results never break the laws of physics. This allows scientists to study complex materials (like new alloys or biological tissues) with much greater confidence and accuracy than before.

1. Problem Statement

The paper addresses the numerical simulation of the Cahn–Hilliard (CH) equations, a fourth-order nonlinear partial differential equation used to model phase separation in binary mixtures. The specific challenges addressed are:

Degenerate Mobility: The mobility function $M(\psi) = \max\{1-\psi^2, 0\}$ vanishes when the phase field $\psi$ reaches $\pm 1$ . This degeneracy complicates the analysis of existence and stability.
Discrete Maximum Principle (DMP): For the quartic double-well potential, the continuous equation does not strictly satisfy a strong maximum principle, but physical solutions are expected to remain bounded within $[-1, 1]$ . Standard Discontinuous Galerkin (DG) and Finite Element Methods (FEM) often violate this bound in numerical simulations, leading to non-physical results.
Arbitrary Polynomial Order: Existing provably bound-preserving schemes are often limited to low-order approximations (e.g., piecewise constant, $p=0$ ) or require complex multi-step procedures. The authors aim to construct a scheme that preserves bounds for arbitrary polynomial order $p$ .
Regularization: Previous approaches often required "regularizing" the mobility (e.g., $M_\delta = \max\{M, \delta\}$ ) to ensure coercivity, which compromises the strict degeneracy of the physical model.

2. Methodology

The authors propose a novel Structure-Preserving Discontinuous Galerkin (DG) scheme called SWIPDP (Symmetric Weighted Interior Penalty Diffusive Projection), combined with a Zhang–Shu scaling limiter (SWIPDP-L).

A. Spatial Discretization (SWIP Formulation)

Cell-Averaged Mobility: A key innovation is evaluating the degenerate mobility $M(\psi)$ $M (ψ)$ at the cell averages ( $\Pi_0 \psi$ $Π_{0} ψ$ ) rather than as a piecewise polynomial.
- This allows the scheme to strictly enforce $M=0$ where the cell average reaches $\pm 1$ , preserving the physical degeneracy without regularization.
- It leads to a coercivity constant that is independent of the mobility value, removing the need for ad-hoc regularization parameters.
Symmetric Weighted Interior Penalty (SWIP): The scheme uses a SWIP formulation for the diffusive flux. The weights are chosen based on the harmonic average of the mobility at cell interfaces. This ensures that if a cell average is at a global maximum/minimum, the flux across that interface is strictly zero, preventing the solution from exceeding bounds.
Time Discretization: The scheme uses a semi-implicit time-stepping approach with Eyre splitting for the nonlinear potential $W(\psi)$ . This splits the potential into convex and concave parts to ensure unconditional energy stability for the unlimited scheme.

B. The Scaling Limiter

To enforce the point-wise bound $|\psi| \le 1$ for $p > 0$ , the authors apply the Zhang–Shu scaling limiter after each time step.

The limiter rescales the polynomial solution on each element to ensure it stays within $[-1, 1]$ while preserving the cell average.
Crucially, the paper proves that the cell averages of the SWIPDP scheme naturally satisfy the bound $|\Pi_0 \psi| \le 1$ (Theorem 4.5). This ensures the limiter has a valid initialization and does not need to alter the cell averages, which is essential for mass conservation.

C. Theoretical Analysis

Existence: The authors prove the existence of discrete solutions for both the unlimited (SWIPDP) and limited (SWIPDP-L) schemes using a Leray–Schauder fixed-point argument. They handle the degenerate limit by passing from a regularized problem ( $\delta > 0$ ) to the degenerate case ( $\delta \to 0$ ).
Coercivity: A rigorous proof establishes that the bilinear form is coercive over a specific subspace, with constants independent of the mobility.
Energy Stability:
- For the unlimited scheme, monotone energy dissipation is proven for all $p$ .
- For the limited scheme ( $p > 0$ ), a formal proof of monotone energy dissipation is not established due to the interaction between the limiter and the energy functional. However, the scheme is shown to be energy-bounded, and numerical experiments confirm monotone dissipation in practice.
Mass Conservation: Mass is conserved exactly for the unlimited scheme. For the limited scheme, mass is conserved because the limiter preserves cell averages (Lemma 4.4).

3. Key Contributions

Arbitrary Order DMP: The first DG scheme to provide a provable discrete maximum principle for the Cahn–Hilliard equation with degenerate mobility at arbitrary polynomial order $p$ .
Strict Degeneracy without Regularization: By evaluating mobility at cell averages, the scheme avoids the need for regularization parameters ( $\delta$ ), which are often required in previous works (e.g., [25, 26]) to ensure coercivity.
Simplified Limiter Workflow: Unlike other high-order bound-preserving schemes (e.g., [27]) that require solving constrained convex optimization problems or multi-step procedures, the proposed scheme requires only a single application of the Zhang–Shu limiter after the time step.
Existence Proof: A rigorous proof of the existence of solutions for the degenerate CH equations using a DG framework, overcoming the lack of coercivity in the degenerate regions.

4. Numerical Results

The authors validate the scheme using the Dune-Fem software package with several test cases:

Convergence Rates: Numerical experiments confirm optimal convergence rates of order $p+1$ in the $L^2$ norm and order $p$ in the $H^1$ norm for both $p=0$ and $p=1$ .
Spinodal Decomposition: Simulations demonstrate that the scheme correctly captures phase separation patterns.
Bound Preservation:
- For $p=0$ , the scheme strictly satisfies the maximum principle (bounds are never violated).
- For $p=1$ , the scheme satisfies the principle up to machine precision (floating-point errors). The deviation scales linearly with the solver tolerance, confirming the theoretical robustness.
Energy and Mass:
- Mass is conserved to machine precision.
- Energy dissipation is observed to be monotone in all numerical experiments, even for $p=1$ where a formal proof is lacking.
Merging Bubbles: The scheme successfully simulates the merging of two bubbles without oscillations or unphysical overshoots, outperforming previous SWIP-L schemes which showed oscillations in similar setups.

5. Significance and Outlook

This work represents a significant advancement in the numerical analysis of phase-field models. By combining a structure-preserving DG formulation with a scaling limiter, the authors have created a robust framework that:

Eliminates the trade-off between high-order accuracy and physical bound preservation.
Removes the need for regularization, ensuring the numerical model adheres strictly to the physical degeneracy of the mobility.
Provides a foundation for more complex coupled systems (e.g., Cahn–Hilliard–Navier–Stokes) where maintaining physical bounds is critical for stability.

The primary limitation noted is the lack of a formal proof for monotone energy dissipation in the limited scheme for $p > 0$ , though numerical evidence strongly supports it. Future work aims to close this theoretical gap and extend the method to adaptive $hp$-meshes and coupled fluid-structure problems.