Mesh: Operators: Cahn Hilliard¶

The “Cahn-Hilliard” equation separates a field $$\phi$$ into 0 and 1 with smooth transitions.

$\frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \left( \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi \right)$

Where $$f$$ is the energy function $$f = ( a^2 / 2 )\phi^2(1 - \phi)^2$$ which drives $$\phi$$ towards either 0 or 1, this competes with the term $$\epsilon^2 \nabla^2 \phi$$ which is a diffusion term that creates smooth changes in $$\phi$$. The equation can be factored:

$\begin{split}\frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \psi \\ \psi = \frac{\partial^2 f}{\partial \phi^2} (\phi - \phi^{\text{old}}) + \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi\end{split}$

Here we will need the derivatives of $$f$$:

$\frac{\partial f}{\partial \phi} = (a^2/2)2\phi(1-\phi)(1-2\phi) \frac{\partial^2 f}{\partial \phi^2} = (a^2/2)2[1-6\phi(1-\phi)]$

The implementation below uses backwards Euler in time with an exponentially increasing time step. The initial $$\phi$$ is a normally distributed field with a standard deviation of 0.1 and mean of 0.5. The grid is 60x60 and takes a few seconds to solve ~130 times. The results are seen below, and you can see the field separating as the time increases. Out:

0 0.006737946999085467
10 0.09636267449939614
20 0.24412886910986079
30 0.4877541372545481
40 0.8894242989247158
50 1.5516664382758794
60 2.643519139778099
70 4.443679913216204
80 7.411643271063606
90 12.304987589805194
100 20.37274845297408
110 33.67423739500265
120 55.604685145707705


from __future__ import print_function
from SimPEG import Mesh, Utils, Solver
import numpy as np
import matplotlib.pyplot as plt

def run(plotIt=True, n=60):

np.random.seed(5)

# Here we are going to rearrange the equations:

# (phi_ - phi)/dt = A*(d2fdphi2*(phi_ - phi) + dfdphi - L*phi_)
# (phi_ - phi)/dt = A*(d2fdphi2*phi_ - d2fdphi2*phi + dfdphi - L*phi_)
# (phi_ - phi)/dt = A*d2fdphi2*phi_ + A*( - d2fdphi2*phi + dfdphi - L*phi_)
# phi_ - phi = dt*A*d2fdphi2*phi_ + dt*A*(- d2fdphi2*phi + dfdphi - L*phi_)
# phi_ - dt*A*d2fdphi2 * phi_ =  dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
# (I - dt*A*d2fdphi2) * phi_ =  dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
# (I - dt*A*d2fdphi2) * phi_ =  dt*A*dfdphi - dt*A*d2fdphi2*phi - dt*A*L*phi_ + phi
# (dt*A*d2fdphi2 - I) * phi_ =  dt*A*d2fdphi2*phi + dt*A*L*phi_ - phi - dt*A*dfdphi
# (dt*A*d2fdphi2 - I - dt*A*L) * phi_ =  (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi

h = [(0.25, n)]
M = Mesh.TensorMesh([h, h])

# Constants
D = a = epsilon = 1.
I = Utils.speye(M.nC)

# Operators
A = D * M.faceDiv * M.cellGrad
L = epsilon**2 * M.faceDiv * M.cellGrad

duration = 75
elapsed = 0.
dexp = -5
phi = np.random.normal(loc=0.5, scale=0.01, size=M.nC)
ii, jj = 0, 0
PHIS = []
capture = np.logspace(-1, np.log10(duration), 8)
while elapsed < duration:
dt = min(100, np.exp(dexp))
elapsed += dt
dexp += 0.05

dfdphi = a**2 * 2 * phi * (1 - phi) * (1 - 2 * phi)
d2fdphi2 = Utils.sdiag(a**2 * 2 * (1 - 6 * phi * (1 - phi)))

MAT = (dt*A*d2fdphi2 - I - dt*A*L)
rhs = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi
phi = Solver(MAT)*rhs

if elapsed > capture[jj]:
PHIS += [(elapsed, phi.copy())]
jj += 1
if ii % 10 == 0:
print(ii, elapsed)
ii += 1

if plotIt:
fig, axes = plt.subplots(2, 4, figsize=(14, 6))
axes = np.array(axes).flatten().tolist()
for ii, ax in zip(np.linspace(0, len(PHIS)-1, len(axes)), axes):
ii = int(ii)
M.plotImage(PHIS[ii], ax=ax)
ax.axis('off')
ax.set_title('Elapsed Time: {0:4.1f}'.format(PHIS[ii]))

if __name__ == '__main__':
run()
plt.show()


Total running time of the script: ( 0 minutes 20.724 seconds)

Gallery generated by Sphinx-Gallery