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