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
                idx = idx%N
    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]$$

  • $\begingroup$ 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? $\endgroup$ – stefano Aug 10 '20 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.