Imposing no flux boundary conditions on variants of the Cahn-Hilliard equation using Finite Differences in Python

Question

I have been looking into simulations of phase separation in variants of the Cahn-Hilliard system and have been running into issues with implementing no flux boundary conditions on certain variants.

The Cahn-Hilliard (C-H) Equation here is: $$\partial_t \phi = -\nabla\cdot\textbf{J}$$ $$\textbf{J} = -\nabla\mu$$ $$\mu = A(\phi^3-\phi) - K\nabla^2\phi$$ where $\phi$ is a nondimensional density with $\phi=1$ representing a high density reference phase and $\phi=-1$ representing a low density reference phase. $A$ and $K$ are constants, $\mu$ is the chemical potential, and $\textbf{J}$ is the total flux.

The boundary conditions are: $$\nabla\phi = 0$$ $$\textbf{J} = -\nabla\mu = 0$$ The first being a Neumann condition representing no diffusive flux into the boundary and the second being a Robin condition representing no total flux into the boundary.

I have run finite difference simulations of C-H with these conditions just fine, but if I introduce certain contributions to the total flux, the system no longer conserves mass when a boundary is introduced (though it does given periodic boundaries). An example is: $$\partial_t \phi = -\nabla\cdot\textbf{J}$$ $$\textbf{J} = -\nabla\mu + D\nabla\phi$$ $$\mu = A(\phi^3-\phi) - K\nabla^2\phi$$ Again, the boundary conditions are: $$\nabla\phi = 0$$ $$\textbf{J} = -\nabla\mu + D\nabla\phi= 0$$ Given the first (Neumann) boundary condition must hold, the second (Robin) boundary condition should become $$\nabla\mu = 0$$ Though I am fairly certain I have implemented this the same way as in the case of just the C-H system, putting in a boundary appears to result in a gradual increase in total mass.

My implementation in python is outlined below as well as the construction of a matrix operator for a finite difference laplacian. This is given a 3-point stencil such that, given lattice spacing $\Delta$ in 1D: $$\nabla^2\phi_i = \frac{\phi_{i+1}-2\phi_i+\phi_{i-1}}{\Delta^2}$$ And given the $\nabla\phi=0$ boundary condition, using a central finite difference scheme we can set: $$\phi_{-1}=\phi_{1}$$ $$\phi_{N} = \phi_{N-2}$$ Since $\mu$ is subject to an analogous boundary condition, I can apply the same operator to $\mu$.

import numpy as np
from scipy.sparse import csr_matrix

dt = 0.02
h = int(10000/dt)
N = 256

def make_laplacian(N, bounds):
    ind = []
    dat = []
    indptr = [0]
    stencil = (1,-2,1)
    for i in range(N):
        for j in range(3):
            idx = i+j-1  #span i-1, i, i+1
            if bounds == True:
                if idx == -1:
                    idx = 1
                elif idx == N:
                    idx = N-2
            else:
                idx = idx%N
            ind.append(idx)
            dat.append(stencil[j])
        indptr.append(len(ind))
    return csr_matrix((dat,ind,indptr),shape = (N,N))

phi = 0.2*(0.5 - np.random.random(N))
D2 = make_laplacian(N,True)
for i in range(h):
    u = phi**3-phi - D2@phi
    phi += dt*(D2@u)-dt*(D2@phi)
    if i%5000==0:
        print(f'Avg phi: {np.mean(phi)}')

Answer 1

The short answer is that you need $$\phi_{-1} = \phi_0$$ $$\phi_N = \phi_{N-1}$$ to impose $\nabla\phi=0$.

A quick check by making the following change

if idx == -1:
    idx = 0
elif idx == N:
    idx = N-1

in the code, you have posted shows that the average $\phi$ remains constant up to 14 decimal places.

To see why this is the correct boundary condition even when you are using central difference, consider the stencil for $\nabla^2\phi$ in 1D:

$$\nabla^2\phi_i = \frac{\phi_{i+1} - 2\phi_i + \phi_{i-1}}{\Delta^2}$$

In here, we are secretly using central difference on a grid spacing of $\Delta/2$ to calculate the derivative of $\nabla \phi$ (see here):

$$\nabla^2\phi_i = \frac{1}{\Delta} \left(\frac{\phi_{i+1} - \phi_i}{\Delta} - \frac{\phi_i - \phi_{i-1}}{\Delta}\right) = \frac{1}{2 (\Delta/2)}(\nabla\phi_{i+1/2} - \nabla\phi_{i-1/2})$$

(If we were using the central difference derivatives on a grid spacing of $\Delta$, then this second derivative would turn out to be $(\phi_{i+2} - 2\phi_i + \phi_{i-2})/(4\Delta^2)$, which would be less accurate.)

So, the derivatives that you need to set to zero for Neumann boundary conditions are $(\phi_N-\phi_{N-1})/\Delta$ and $(\phi_0-\phi_{-1})/\Delta$.

An easy way to check what's going wrong with mass conservation in finite difference schemes is to explicitly compute the mass and check the residue:

$$\phi^{t+1}_i = \phi^t_i + \frac{\mathrm{d}t}{\Delta^2}(\mu_{i+1} - 2\mu_i + \mu_{i-1})$$

$$\implies \sum_{i=0}^{N-1}\phi^{t+1}_i = \sum_{i=0}^{N-1}\phi^t_i + \frac{\mathrm{d}t}{\Delta^2}\sum_{i=0}^{N-1}(\mu_{i+1} - 2\mu_i + \mu_{i-1})$$

$$\implies M^{t+1} = M^t + \frac{\mathrm{d}t}{\Delta^2}\left[(\mu_{N} - \mu_{N-1}) - (\mu_0- \mu_{-1})\right]$$