# 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.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 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): """ Mesh: Operators: Cahn Hilliard ============================== This example is based on the example in the FiPy_ library. Please see their documentation for more information about the Cahn-Hilliard equation. The "Cahn-Hilliard" equation separates a field \\\$$\\\\phi \\\$$ into 0 and 1 with smooth transitions. .. math:: \\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: .. math:: \\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 Here we will need the derivatives of \\\$$f \\\$$: .. math:: \\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. .. _FiPy: http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html """ 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][1], ax=ax) ax.axis('off') ax.set_title('Elapsed Time: {0:4.1f}'.format(PHIS[ii][0])) if __name__ == '__main__': run() plt.show()