From b2aab4163b655b93c1465012b9e5e2024e308e05 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 26 Nov 2015 13:56:04 -0700 Subject: [PATCH] New examples for Mesh --- .../Examples/Mesh_Operators_CahnHilliard.py | 105 ++++++++++++++++++ SimPEG/Examples/Mesh_QuadTree_FaceDiv.py | 49 ++++++++ docs/examples/Mesh_Operators_CahnHilliard.rst | 57 ++++++++++ docs/examples/Mesh_QuadTree_FaceDiv.rst | 26 +++++ 4 files changed, 237 insertions(+) create mode 100644 SimPEG/Examples/Mesh_Operators_CahnHilliard.py create mode 100644 SimPEG/Examples/Mesh_QuadTree_FaceDiv.py create mode 100644 docs/examples/Mesh_Operators_CahnHilliard.rst create mode 100644 docs/examples/Mesh_QuadTree_FaceDiv.rst diff --git a/SimPEG/Examples/Mesh_Operators_CahnHilliard.py b/SimPEG/Examples/Mesh_Operators_CahnHilliard.py new file mode 100644 index 00000000..8e42e617 --- /dev/null +++ b/SimPEG/Examples/Mesh_Operators_CahnHilliard.py @@ -0,0 +1,105 @@ +from SimPEG import * + +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: + import matplotlib.pyplot as plt + 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) + out = M.plotImage(PHIS[ii][1],ax=ax) + ax.axis('off') + ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0]) + plt.show() + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py b/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py new file mode 100644 index 00000000..5bd67929 --- /dev/null +++ b/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py @@ -0,0 +1,49 @@ +from SimPEG import * + +def run(plotIt=True, n=60): + """ + Mesh: QuadTree: FaceDiv + ======================= + + + + """ + + + M = Mesh.TreeMesh([[(1,16)],[(1,16)]], levels=4) + M._refineCell([0,0,0]) + M._refineCell([0,0,1]) + M._refineCell([4,4,2]) + M.__dirty__ = True + M.number() + + + if plotIt: + import matplotlib.pyplot as plt + fig, axes = plt.subplots(2,1,figsize=(10,10)) + + M.plotGrid(cells=True, nodes=False, ax=axes[0]) + axes[0].axis('off') + axes[0].set_title('Simple QuadTree Mesh') + axes[0].set_xlim([-1,17]) + axes[0].set_ylim([-1,17]) + + for ii, loc in zip(range(M.nC),M.gridCC): + axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r') + + axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>') + for ii, loc in zip(range(M.nFx),M.gridFx): + axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g') + + axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^') + for ii, loc in zip(range(M.nFy),M.gridFy): + axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m') + + axes[1].spy(M.faceDiv) + axes[1].set_title('Face Divergence') + axes[1].set_ylabel('Cell Number') + axes[1].set_xlabel('Face Number') + plt.show() + +if __name__ == '__main__': + run() diff --git a/docs/examples/Mesh_Operators_CahnHilliard.rst b/docs/examples/Mesh_Operators_CahnHilliard.rst new file mode 100644 index 00000000..9786e911 --- /dev/null +++ b/docs/examples/Mesh_Operators_CahnHilliard.rst @@ -0,0 +1,57 @@ +.. _examples_Mesh_Operators_CahnHilliard: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +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 + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_Operators_CahnHilliard.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_Operators_CahnHilliard.py + :language: python + :linenos: diff --git a/docs/examples/Mesh_QuadTree_FaceDiv.rst b/docs/examples/Mesh_QuadTree_FaceDiv.rst new file mode 100644 index 00000000..6bfdd47f --- /dev/null +++ b/docs/examples/Mesh_QuadTree_FaceDiv.rst @@ -0,0 +1,26 @@ +.. _examples_Mesh_QuadTree_FaceDiv: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Mesh: QuadTree: FaceDiv +======================= + + + + + +.. plot:: + + from SimPEG import Examples + Examples.Mesh_QuadTree_FaceDiv.run() + +.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_FaceDiv.py + :language: python + :linenos: