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

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 = 
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)}')


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]$$

• I see how it checks out from the residual, and I verified that the code works with your modification. However, it's common in FD schemes to use a ghost point to implement Neumann BC to second order accuracy, which is what Aaron seems to be doing. You can find it described everywhere, for instance on page 4 of math.uci.edu/~chenlong/226/FDM.pdf Do you have an idea as to why it doesn't work? Is it being done incorrectly here? Aug 10 '20 at 15:17