Task 4: A first IMEX solver for the Cahn-Hilliard equation

Now we have finally arrived at the point where we can start implementing a first solver for the Cahn-Hilliard equation. You are tasked with extending (parts of) your transient biharmonic solver from the previous task to include the additional non-linearity \(\Delta f(u)\) arising from the Cahn-Hilliard equation. The non-linearity poses several challenge for the time-stepping scheme:

• As an explicit time-stepping scheme needs to satisfy severe time-step restrictions when combined with the transient biharmonic operator (see previous task), an implicit time-stepping scheme seems to be more appropriate. On the other hand, the non-linearity \(\Delta f(u)\) is computationally expensive to solve implicitly, as it would require the solution of a non-linear system of equations at each time step.

• Moreover, as the product of functions translates into convolution in the Fourier space, the non-linearity \( f(u)\) translates into a convolution of the form \(\widehat{u^3}(k) - \widehat{u}(k) = \widehat{u}\ast\widehat{u}\ast \widehat{u}(k) - \widehat{u}(k)\), which is extremely expensive to compute!

To address these challenge, we need to briefly discuss some import concepts when solving the Cahn-Hilliard equation.

IMEX time-stepping scheme based on the Implicit/Explicit Euler method

The first challenge can be addressed by using an Implicit-Explicit (IMEX) time-stepping scheme. The IMEX approach is best explained by looking at a general ODE system of the form

\[ U_t + \mathbf{L}U = \mathbf{N}(U), \quad U(0) = U_0, \]

Again, we start from a equidistant time grid \(t_n = n\tau\) with \(\tau = (T-t_0)/N_t\) and denote the numerical solution at time \(t_n\) by \(U^n\). Starting from the initial condition \(U^0\) at time \(t_0\), and given the solution \(U^n\) at time \(t_n\), the new solution \(U^{n+1}\) at time step \(t_{n+1} = t_n + \tau\) is then computed as

\[ U^{n+1} + \tau \mathbf{L}U^{n+1} = U^n + \tau \mathbf{N}(U^{n}). \]

This explains the name “IMEX” (Implicit-Explicit) time-stepping scheme, as the linear parts of the equation are treated implicitly, while the non-linear part is treated explicitly. In particular, for \(\mathbf{N} = 0\), the scheme reduces to the Implicit Euler method, while for \(\mathbf{L} = 0\), the scheme reduces to the Explicit Euler method.

Note that this idea will be used to solve the Cahn-Hilliard equation in the Fourier space, so below you have to think of \(U^n\) being \(\widehat{u}^n\) and \(\mathbf{L}\) being the discrete Fourier transform of the Biharmonic operator (potentially plus some lower terms).

Pseudo-spectral methods for nonlinear PDEs

To address the second challenge, we will implement a so-called pseudo-spectral method. The idea is to compute the non-linear term in the physical space, while the linear terms are computed in the Fourier space. This allows for an efficient computation of the non-linear term, while still benefiting from the high accuracy of the Fourier space method for the linear terms. So to compute \(\widehat{N(u^n)}\) at \(t_{n+1}\), we start from the (DFT of the) solution of the previous time step \(\widehat{u}^n\) and proceed as follows:

\[ \widehat{u}^n \overset{\mathcal{F}^{-1}}{\mapsto} u^n \mapsto N(u^n) \overset{\mathcal{F}}{\mapsto} \widehat{N(u^n)}. \]

So in contrast to the previous task, where we only transformed the solution back to the physical space for “post-processing” tasks such as visualization or error computations, we now need to transform the solution back and forth for the actual solution computation. Hence, the name pseudo-spectral method is frequently used to refer to such kind of methods, as we cannot stay in the Fourier space for the entire computation.

Convex splitting of the Cahn-Hilliard equation

Finally, we need to specify how exactly we want to split the linear parts from the non-linear parts of the Cahn-Hilliard equation. Note that the the term

\[ \Delta f(u) = \Delta (u^3 -u) \]

in the Cahn-Hilliard equation actually contains a linear part \(-\Delta u\) and a non-linear part \(\Delta u^3\).

If we decide to treat all linear terms implicitly and all non-linear terms explicitly, we can split the Cahn-Hilliard equation as follows:

(85)\[\begin{align} \dfrac{u^{n+1} - u^n}{\tau} + \kappa \Delta^2 u^{n+1} + \Delta u^{n+1} = \Delta (u^n)^3 \end{align}\]

Unfortunately, the \(+\Delta u^{n+1}\) has an unfortunate sign, which we see when transforming this equation to the Fourier space:

(86)\[\begin{align} \widehat{u}^{n+1} + \tau (\kappa \widetilde{\mathbf{k}}^4 \widehat{u}^{n+1} - \widetilde{\mathbf{k}}^2 \widehat{u}^{n+1}) = \widehat{u}^{n} + \tau \widetilde{\mathbf{k}}^2 \widehat{(u^n)^3} \end{align}\]

Consequently, this solution method can get unstable if \(1 + \kappa \widetilde{\mathbf{k}}^4 - \widetilde{\mathbf{k}}^2 < - 1\) which can easily happen for reasonable values of \(\kappa\) and \(\mathbf{k}\).

As a remedy, it is common to split the \(\Delta\) slightly differently, starting from a splitting parameter \(a\) to obtain

\[ \Delta f(u) = \Delta (u^3 -u) = \underbrace{a\Delta u}_{f_1(u)} + \underbrace{\Delta u^3 -(1+a) \Delta u}_{f_2(u)}. \]

where \(f_1(u)\) is treated implicitly and \(f_2(u)\) is treated explicitly. Note that for \(a = -1\), we recover the original complete splitting. But typically, we choose at least \(a \geqslant 0\) to avoid the mentioned sign issue in the Fourier space, but for reasons we don’t have time to discuss here, one chooses \( a \sim 1.5\). For \(a \geqslant 2\), this splitting results from a “convex” splitting of the free energy functional.

The result is then the following IMEX scheme for the Cahn-Hilliard equation:

(87)\[\begin{align} \dfrac{u^{n+1} - u^n}{\tau} +\kappa \Delta^2 u^{n+1} - a \Delta u^{n+1} = \Delta \bigl( (u^n)^3 - (1+a) u^n\bigr) \end{align}\]

which of course needs to be translated to the Fourier space.

Computational tasks

1) Before you start implementing, please provide a brief mathematical description of the resulting numerical scheme, in particular, describe how a new solution is computed in the Fourier space from the previous solution for each time step. To ensure that you later can verify the correctness of your implementation using the method of manufactured solutions, make sure that your solver can solve the Cahn-Hilliard equation with a source term \(g(x,y,t)\), i.e., the equation to be solved should be of the form

(88)\[\begin{alignat}{2} \partial_t u - \nabla\cdot (M \nabla (f(u) - \kappa \Delta u )) &= g & &\quad \text{ on } \Omega \times (0,T) \end{alignat}\]

Treat the source term \(g\) implicitly in the time-stepping scheme (similar to the linear terms).

2) Implement the resulting IMEX scheme for the Cahn-Hilliard equation. For the implementation, the solver interface should meet similar specifications as in the previous task, i.e., the solvers should be implemented as a generator function using the yield statement to return the discrete fourier transform of the solution at each time step, and the current time. The generator function should have the following signature:

def cahn_hilliard_backward_euler(*, 
                                kappa, 
                                X, Y, U0, 
                                t0, T, Nt,
                                g, 
                                alpha=1.5):
    """
    Implements the Cahn-Hilliard equation solver using the backward Euler method 
    with a convex-concave splitting approach.

    Parameters:
    -----------
    kappa : float
        Diffusion coefficient for the biharmonic operator.
    X : ndarray
        2D array representing the x-coordinates of the grid.
    Y : ndarray
        2D array representing the y-coordinates of the grid.
    U0 : ndarray
        Initial condition for the solution.
    t0 : float
        Initial time.
    T : float
        Final time.
    Nt : int
        Number of time steps.
    g : callable or None
        Source term as a function of (X, Y, t). If None, no source term is applied.
    alpha : float, optional
        Convex-concave splitting parameter. Default is 1.5.

    Yields:
    -------
    tuple: A tuple containing the discrete Fourier transform of U at t, and the current time t.

    """

3) As before, we want you to verify the correctness of your implementation by comparing the numerical solution to the exact solution for a simple test case. To this end, we ask you to run a convergence study similar to the one in the previous task.

To manufacture a solution, start from the exact solution

\[ u_{\mathrm{ex}} = \sin(x)\cos(y)\exp(-4 \kappa t) \]

Clearly, this solution will not satisfy the Cahn-Hilliard equation with a non-zero source term due to the additional non-linearity. Instead, compute the corresponding right-hand side \(g\) for the Cahn-Hilliard equation. This might be a good time to automate the computation of the right-hand side using symbolic computation tools such as sympy (see lecture notes)!

Now for \(\Omega = [0,16\pi)^2\), \(N_x = N_y = 64\), \(t_0 = 0\), \(T = 1\) and \(\kappa \in \{1.0, 0.01\}\), run convergence studies for the IMEX solver using \(N_t \in \{100, 200, 400, 800, 1600, 3200\}\), where you tabulate the \(L^{\infty}L^{\infty}\) error against the number of number of time steps \(N_t\). Discuss briefly your results.