examples.cahnHilliard package¶

Submodules¶

examples.cahnHilliard.mesh2D module¶

The spinodal decomposition phenomenon is a spontaneous separation of an initially homogeneous mixture into two distinct regions of different properties (spin-up/spin-down, component A/component B). It is a “barrierless” phase separation process, such that under the right thermodynamic conditions, any fluctuation, no matter how small, will tend to grow. This is in contrast to nucleation, where a fluctuation must exceed some critical magnitude before it will survive and grow. Spinodal decomposition can be described by the “Cahn-Hilliard” equation (also known as “conserved Ginsberg-Landau” or “model B” of Hohenberg & Halperin)

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

where $\phi$ is a conserved order parameter, possibly representing alloy composition or spin. The double-well free energy function $f = (a^2/2) \phi^2 (1 - \phi)^2$ penalizes states with intermediate values of $\phi$ between 0 and 1. The gradient energy term $\epsilon^2 \nabla^2\phi$ , on the other hand, penalizes sharp changes of $\phi$ . These two competing effects result in the segregation of $\phi$ into domains of 0 and 1, separated by abrupt, but smooth, transitions. The parameters $a$ and $\epsilon$ determine the relative weighting of the two effects and $D$ is a rate constant.

We can simulate this process in FiPy with a simple script:

>>> from fipy import CellVariable, Grid2D, GaussianNoiseVariable, TransientTerm, DiffusionTerm, ImplicitSourceTerm, LinearLUSolver, Viewer
>>> from fipy.tools import numerix

(Note that all of the functionality of NumPy is imported along with FiPy, although much is augmented for FiPy's needs.)

>>> if __name__ == "__main__":
...     nx = ny = 1000
... else:
...     nx = ny = 20
>>> mesh = Grid2D(nx=nx, ny=ny, dx=0.25, dy=0.25)
>>> phi = CellVariable(name=r"$\phi$", mesh=mesh)

We start the problem with random fluctuations about $\phi = 1/2$

>>> phi.setValue(GaussianNoiseVariable(mesh=mesh,
...                                    mean=0.5,
...                                    variance=0.01))

FiPy doesn’t plot or output anything unless you tell it to:

>>> if __name__ == "__main__":
...     viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)

For FiPy, we need to perform the partial derivative $\partial f/\partial \phi$ manually and then put the equation in the canonical form by decomposing the spatial derivatives so that each Term is of a single, even order:

$\frac{\partial \phi}{\partial t} = \nabla\cdot D a^2 \left[ 1 - 6 \phi \left(1 - \phi\right)\right] \nabla \phi- \nabla\cdot D \nabla \epsilon^2 \nabla^2 \phi.$

FiPy would automatically interpolate D * a**2 * (1 - 6 * phi * (1 - phi)) onto the faces, where the diffusive flux is calculated, but we obtain somewhat more accurate results by performing a linear interpolation from phi at cell centers to PHI at face centers. Some problems benefit from non-linear interpolations, such as harmonic or geometric means, and FiPy makes it easy to obtain these, too.

>>> PHI = phi.arithmeticFaceValue
>>> D = a = epsilon = 1.
>>> eq = (TransientTerm()
...       == DiffusionTerm(coeff=D * a**2 * (1 - 6 * PHI * (1 - PHI)))
...       - DiffusionTerm(coeff=(D, epsilon**2)))

Because the evolution of a spinodal microstructure slows with time, we use exponentially increasing time steps to keep the simulation “interesting”. The FiPy user always has direct control over the evolution of their problem.

>>> dexp = -5
>>> elapsed = 0.
>>> if __name__ == "__main__":
...     duration = 1000.
... else:
...     duration = 1000.

>>> while elapsed < duration:
...     dt = min(100, numerix.exp(dexp))
...     elapsed += dt
...     dexp += 0.01
...     eq.solve(phi, dt=dt, solver=LinearLUSolver())
...     if __name__ == "__main__":
...         viewer.plot()
...     elif (max(phi.globalValue) > 0.7) and (min(phi.globalValue) < 0.3) and elapsed > 10.:
...         break

>>> print((max(phi.globalValue) > 0.7) and (min(phi.globalValue) < 0.3))
True

evolution of Cahn-Hilliard phase separation at t = 30, 100 and 1000

examples.cahnHilliard.mesh2DCoupled module¶

Solve the Cahn-Hilliard problem in two dimensions.

Warning

This formulation has serious performance problems and is not automatically tested. Specifically, for non-trivial mesh sizes, PySparse requires enormous amounts of memory, Trilinos cannot solve the coupled form, and PETSc cannot solve the vector form.

The spinodal decomposition phenomenon is a spontaneous separation of an initially homogeneous mixture into two distinct regions of different properties (spin-up/spin-down, component A/component B). It is a “barrierless” phase separation process, such that under the right thermodynamic conditions, any fluctuation, no matter how small, will tend to grow. This is in contrast to nucleation, where a fluctuation must exceed some critical magnitude before it will survive and grow. Spinodal decomposition can be described by the “Cahn-Hilliard” equation (also known as “conserved Ginsberg-Landau” or “model B” of Hohenberg & Halperin)

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

where $\phi$ is a conserved order parameter, possibly representing alloy composition or spin. The double-well free energy function $f = (a^2/2) \phi^2 (1 - \phi)^2$ penalizes states with intermediate values of $\phi$ between 0 and 1. The gradient energy term $\epsilon^2 \nabla^2\phi$ , on the other hand, penalizes sharp changes of $\phi$ . These two competing effects result in the segregation of $\phi$ into domains of 0 and 1, separated by abrupt, but smooth, transitions. The parameters $a$ and $\epsilon$ determine the relative weighting of the two effects and $D$ is a rate constant.

We can simulate this process in FiPy with a simple script:

>>> from fipy import CellVariable, Grid2D, GaussianNoiseVariable, DiffusionTerm, TransientTerm, ImplicitSourceTerm, Viewer
>>> from fipy.tools import numerix

(Note that all of the functionality of NumPy is imported along with FiPy, although much is augmented for FiPy's needs.)

>>> if __name__ == "__main__":
...     nx = ny = 20
... else:
...     nx = ny = 10
>>> mesh = Grid2D(nx=nx, ny=ny, dx=0.25, dy=0.25)
>>> phi = CellVariable(name=r"$\phi$", mesh=mesh)
>>> psi = CellVariable(name=r"$\psi$", mesh=mesh)

We start the problem with random fluctuations about $\phi = 1/2$

>>> noise = GaussianNoiseVariable(mesh=mesh,
...                               mean=0.5,
...                               variance=0.01).value

>>> phi[:] = noise

FiPy doesn’t plot or output anything unless you tell it to:

>>> if __name__ == "__main__":
...     viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)

We factor the Cahn-Hilliard equation into two 2nd-order PDEs and place them in canonical form for FiPy to solve them as a coupled set of equations.

$\frac{\partial \phi}{\partial t} &= \nabla\cdot D \nabla \psi \\ \psi &= \left(\frac{\partial f}{\partial \phi} - \frac{\partial^2 f}{\partial \phi^2}\phi\right)_{\text{old}} + \frac{\partial^2 f}{\partial \phi^2}\phi - \epsilon^2 \nabla^2 \phi$

The source term in $\psi$ , $\frac{\partial f}{\partial \phi}$ , is expressed in linearized form after Taylor expansion at $\phi=\phi_{\text{old}}$ , for the same reasons discussed in examples.phase.simple. We need to perform the partial derivatives

$\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 \left[1 - 6 \phi(1 - \phi)\right]$

manually.

>>> D = a = epsilon = 1.
>>> dfdphi = a**2 * phi * (1 - phi) * (1 - 2 * phi)
>>> dfdphi_ = a**2 * (1 - phi) * (1 - 2 * phi)
>>> d2fdphi2 = a**2 * (1 - 6 * phi * (1 - phi))
>>> eq1 = (TransientTerm(var=phi) == DiffusionTerm(coeff=D, var=psi))
>>> eq2 = (ImplicitSourceTerm(coeff=1., var=psi)
...        == ImplicitSourceTerm(coeff=d2fdphi2, var=phi) - d2fdphi2 * phi + dfdphi
...        - DiffusionTerm(coeff=epsilon**2, var=phi))
>>> eq3 = (ImplicitSourceTerm(coeff=1., var=psi)
...        == ImplicitSourceTerm(coeff=dfdphi_, var=phi)
...        - DiffusionTerm(coeff=epsilon**2, var=phi))

>>> eq = eq1 & eq2

Because the evolution of a spinodal microstructure slows with time, we use exponentially increasing time steps to keep the simulation “interesting”. The FiPy user always has direct control over the evolution of their problem.

>>> dexp = -5
>>> elapsed = 0.
>>> if __name__ == "__main__":
...     duration = 100.
... else:
...     duration = .5e-1

>>> while elapsed < duration:
...     dt = min(100, numerix.exp(dexp))
...     elapsed += dt
...     dexp += 0.01
...     eq.solve(dt=dt)
...     if __name__ == "__main__":
...         viewer.plot()

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Coupled equations. Press <return> to proceed...")

documentation/generated/examples/mesh2DCoupled.*

These equations can also be solved in FiPy using a vector equation. The variables $\phi$ and $\psi$ are now stored in a single variable

>>> var = CellVariable(mesh=mesh, elementshape=(2,))
>>> var[0] = noise

>>> if __name__ == "__main__":
...     viewer = Viewer(name=r"$\phi$", vars=var[0,], datamin=0., datamax=1.)

>>> D = a = epsilon = 1.
>>> v0 = var[0]
>>> dfdphi = a**2 * v0 * (1 - v0) * (1 - 2 * v0)
>>> dfdphi_ = a**2 * (1 - v0) * (1 - 2 * v0)
>>> d2fdphi2 = a**2 * (1 - 6 * v0 * (1 - v0))

The source terms have to be shaped correctly for a vector. The implicit source coefficient has to have a shape of (2, 2) while the explicit source has a shape (2,)

>>> source = (- d2fdphi2 * v0 + dfdphi) * (0, 1)
>>> impCoeff = d2fdphi2 * ((0, 0),
...                        (1., 0)) + ((0, 0),
...                                    (0, -1.))

This is the same equation as the previous definition of eq, but now in a vector format.

>>> eq = (TransientTerm(((1., 0.),
...                     (0., 0.))) == DiffusionTerm([((0.,          D),
...                                                   (-epsilon**2, 0.))])
...                                   + ImplicitSourceTerm(impCoeff) + source)

>>> dexp = -5
>>> elapsed = 0.

>>> while elapsed < duration:
...     dt = min(100, numerix.exp(dexp))
...     elapsed += dt
...     dexp += 0.01
...     eq.solve(var=var, dt=dt)
...     if __name__ == "__main__":
...         viewer.plot()

>>> print(numerix.allclose(var, (phi, psi)))
True

examples.cahnHilliard.mesh3D module¶

Solves the Cahn-Hilliard problem in a 3D cube

>>> from fipy import CellVariable, Grid3D, Viewer, GaussianNoiseVariable, TransientTerm, DiffusionTerm, DefaultSolver
>>> from fipy.tools import numerix

The only difference from examples.cahnHilliard.mesh2D is the declaration of mesh.

>>> if __name__ == "__main__":
...     nx = ny = nz = 100
... else:
...     nx = ny = nz = 10
>>> mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=0.25, dy=0.25, dz=0.25)
>>> phi = CellVariable(name=r"$\phi$", mesh=mesh)

We start the problem with random fluctuations about $\phi = 1/2$

>>> phi.setValue(GaussianNoiseVariable(mesh=mesh,
...                                    mean=0.5,
...                                    variance=0.01))

FiPy doesn’t plot or output anything unless you tell it to:

>>> if __name__ == "__main__":
...     viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)

For FiPy, we need to perform the partial derivative $\partial f/\partial \phi$ manually and then put the equation in the canonical form by decomposing the spatial derivatives so that each Term is of a single, even order:

$\frac{\partial \phi}{\partial t} = \nabla\cdot D a^2 \left[ 1 - 6 \phi \left(1 - \phi\right)\right] \nabla \phi- \nabla\cdot D \nabla \epsilon^2 \nabla^2 \phi.$

FiPy would automatically interpolate D * a**2 * (1 - 6 * phi * (1 - phi)) onto the Faces, where the diffusive flux is calculated, but we obtain somewhat more accurate results by performing a linear interpolation from phi at Cell centers to PHI at Face centers. Some problems benefit from non-linear interpolations, such as harmonic or geometric means, and FiPy makes it easy to obtain these, too.

>>> PHI = phi.arithmeticFaceValue
>>> D = a = epsilon = 1.
>>> eq = (TransientTerm()
...       == DiffusionTerm(coeff=D * a**2 * (1 - 6 * PHI * (1 - PHI)))
...       - DiffusionTerm(coeff=(D, epsilon**2)))

Because the evolution of a spinodal microstructure slows with time, we use exponentially increasing time steps to keep the simulation “interesting”. The FiPy user always has direct control over the evolution of their problem.

>>> dexp = -5
>>> elapsed = 0.
>>> if __name__ == "__main__":
...     duration = 1000.
... else:
...     duration = 1e-2
>>> while elapsed < duration:
...     dt = min(100, numerix.exp(dexp))
...     elapsed += dt
...     dexp += 0.01
...     eq.solve(phi, dt=dt, solver=DefaultSolver(precon=None))
...     if __name__ == "__main__":
...         viewer.plot()

snapshot of Cahn-Hilliard phase separation in 3D with cutaway

examples.cahnHilliard.sphere module¶

Solves the Cahn-Hilliard problem on the surface of a sphere.

This phenomenon can occur on vesicles (http://www.youtube.com/watch?v=kDsFP67_ZSE).

>>> from fipy import CellVariable, Gmsh2DIn3DSpace, GaussianNoiseVariable, Viewer, TransientTerm, DiffusionTerm, DefaultSolver
>>> from fipy.tools import numerix

The only difference from examples.cahnHilliard.mesh2D is the declaration of mesh.

>>> mesh = Gmsh2DIn3DSpace('''
...     radius = 5.0;
...     cellSize = 0.3;
...
...     // create inner 1/8 shell
...     Point(1) = {0, 0, 0, cellSize};
...     Point(2) = {-radius, 0, 0, cellSize};
...     Point(3) = {0, radius, 0, cellSize};
...     Point(4) = {0, 0, radius, cellSize};
...     Circle(1) = {2, 1, 3};
...     Circle(2) = {4, 1, 2};
...     Circle(3) = {4, 1, 3};
...     Line Loop(1) = {1, -3, 2} ;
...     Ruled Surface(1) = {1};
...
...     // create remaining 7/8 inner shells
...     t1[] = Rotate {{0,0,1},{0,0,0},Pi/2} {Duplicata{Surface{1};}};
...     t2[] = Rotate {{0,0,1},{0,0,0},Pi} {Duplicata{Surface{1};}};
...     t3[] = Rotate {{0,0,1},{0,0,0},Pi*3/2} {Duplicata{Surface{1};}};
...     t4[] = Rotate {{0,1,0},{0,0,0},-Pi/2} {Duplicata{Surface{1};}};
...     t5[] = Rotate {{0,0,1},{0,0,0},Pi/2} {Duplicata{Surface{t4[0]};}};
...     t6[] = Rotate {{0,0,1},{0,0,0},Pi} {Duplicata{Surface{t4[0]};}};
...     t7[] = Rotate {{0,0,1},{0,0,0},Pi*3/2} {Duplicata{Surface{t4[0]};}};
...
...     // create entire inner and outer shell
...     Surface Loop(100)={1,t1[0],t2[0],t3[0],t7[0],t4[0],t5[0],t6[0]};
... ''', overlap=2).extrude(extrudeFunc=lambda r: 1.1 * r) 
>>> phi = CellVariable(name=r"$\phi$", mesh=mesh) 

We start the problem with random fluctuations about $\phi = 1/2$

>>> phi.setValue(GaussianNoiseVariable(mesh=mesh,
...                                    mean=0.5,
...                                    variance=0.01)) 

FiPy doesn’t plot or output anything unless you tell it to: If MayaviClient is available, we can customize the view with a subclass of MayaviDaemon.

>>> if __name__ == "__main__":
...     try:
...         viewer = MayaviClient(vars=phi,
...                               datamin=0., datamax=1.,
...                               daemon_file="examples/cahnHilliard/sphereDaemon.py")
...     except:
...         viewer = Viewer(vars=phi,
...                         datamin=0., datamax=1.,
...                         xmin=-2.5, zmax=2.5)

For FiPy, we need to perform the partial derivative $\partial f/\partial \phi$ manually and then put the equation in the canonical form by decomposing the spatial derivatives so that each Term is of a single, even order:

$\frac{\partial \phi}{\partial t} = \nabla\cdot D a^2 \left[ 1 - 6 \phi \left(1 - \phi\right)\right] \nabla \phi- \nabla\cdot D \nabla \epsilon^2 \nabla^2 \phi.$

FiPy would automatically interpolate D * a**2 * (1 - 6 * phi * (1 - phi)) onto the faces, where the diffusive flux is calculated, but we obtain somewhat more accurate results by performing a linear interpolation from phi at cell centers to PHI at face centers. Some problems benefit from non-linear interpolations, such as harmonic or geometric means, and FiPy makes it easy to obtain these, too.

>>> PHI = phi.arithmeticFaceValue 
>>> D = a = epsilon = 1.
>>> eq = (TransientTerm()
...       == DiffusionTerm(coeff=D * a**2 * (1 - 6 * PHI * (1 - PHI)))
...       - DiffusionTerm(coeff=(D, epsilon**2))) 

Because the evolution of a spinodal microstructure slows with time, we use exponentially increasing time steps to keep the simulation “interesting”. The FiPy user always has direct control over the evolution of their problem.

>>> dexp = -5
>>> elapsed = 0.
>>> if __name__ == "__main__":
...     duration = 1000.
... else:
...     duration = 1e-2
>>> while elapsed < duration:
...     dt = min(100, numerix.exp(dexp))
...     elapsed += dt
...     dexp += 0.01
...     eq.solve(phi, dt=dt, solver=DefaultSolver(precon=None)) 
...     if __name__ == "__main__":
...         viewer.plot()

Cahn-Hilliard phase separation on the surface of a sphere with a rendering of the mesh

examples.cahnHilliard.sphereDaemon module¶

class examples.cahnHilliard.sphereDaemon.SphereDaemon¶

Bases: MayaviDaemon

__base_traits__ = {'application': <traits.ctrait.CTrait object>, 'log_mode': <traits.ctrait.CTrait object>, 'script': <traits.ctrait.CTrait object>, 'start_gui_event_loop': <traits.ctrait.CTrait object>, 'trait_added': <traits.ctrait.CTrait object>, 'trait_modified': <traits.ctrait.CTrait object>}¶

__class_traits__ = {'application': <traits.ctrait.CTrait object>, 'log_mode': <traits.ctrait.CTrait object>, 'script': <traits.ctrait.CTrait object>, 'start_gui_event_loop': <traits.ctrait.CTrait object>, 'trait_added': <traits.ctrait.CTrait object>, 'trait_modified': <traits.ctrait.CTrait object>}¶

__del__()¶

__delattr__(name, /)¶: Implement delattr(self, name).

__getattribute__(name, /)¶: Return getattr(self, name).

__instance_traits__ = {}¶

__listener_traits__ = {'_on_application_gui_started': ('method', {'pattern': 'application.gui:started', 'post_init': False, 'dispatch': 'same'})}¶

__observer_traits__ = {}¶

__prefix_traits__ = {'': <traits.ctrait.CTrait object>, '*': ['_traits_cache_', ''], '_traits_cache_': <traits.ctrait.CTrait object>}¶

__setattr__(name, value, /)¶: Implement setattr(self, name, value).

__view_traits__ = {}¶

classmethod add_class_trait(name, *trait)¶

Adds a named trait attribute to this class.

Also adds the same attribute to all subclasses.

Parameters:

name (str) – Name of the attribute to add.
*trait – A trait or a value that can be converted to a trait using Trait() Trait definition of the attribute. It can be a single value or a list equivalent to an argument list for the Trait() function.

classmethod class_default_traits_view()¶: Returns the name of the default traits view for the class.

classmethod class_editable_traits()¶: Returns an alphabetically sorted list of the names of non-event trait attributes associated with the current class.

classmethod class_trait_names(**metadata)¶

Returns a list of the names of all trait attributes whose definitions match the set of metadata criteria specified.

This method is similar to the traits() method, but returns only the names of the matching trait attributes, not the trait definitions.

Parameters:: **metadata – Criteria for selecting trait attributes.

classmethod class_trait_view(name=None, view_element=None)¶

classmethod class_trait_view_elements()¶

Returns the ViewElements object associated with the class.

The returned object can be used to access all the view elements associated with the class.

classmethod class_traits(**metadata)¶

Returns a dictionary containing the definitions of all of the trait attributes of the class that match the set of metadata criteria.

The keys of the returned dictionary are the trait attribute names, and the values are their corresponding trait definition objects.

If no metadata information is specified, then all explicitly defined trait attributes defined for the class are returned.

Otherwise, the metadata keyword dictionary is assumed to define a set of search criteria for selecting trait attributes of interest. The metadata dictionary keys correspond to the names of trait metadata attributes to examine, and the values correspond to the values the metadata attribute must have in order to be included in the search results.

The metadata values either may be simple Python values like strings or integers, or may be lambda expressions or functions that return True if the trait attribute is to be included in the result. A lambda expression or function must receive a single argument, which is the value of the trait metadata attribute being tested. If more than one metadata keyword is specified, a trait attribute must match the metadata values of all keywords to be included in the result.

Parameters:: **metadata – Criteria for selecting trait attributes.

classmethod class_visible_traits()¶: Returns an alphabetically sorted list of the names of non-event trait attributes associated with the current class, that should be GUI visible

clip_data(src)¶

parse_command_line(argv)¶

Parse command line options.

Parameters:: argv (list of str) – The command line arguments

poll_file()¶

run()¶

This function is called after the GUI has started. Override this to do whatever you want to do as a MayaVi script. If this is not overridden then an empty MayaVi application will be started.

Make sure all other MayaVi specific imports are made here! If you import MayaVi related code earlier you will run into difficulties. Use ‘self.script’ to script the mayavi engine.

classmethod set_trait_dispatch_handler(name, klass, override=False)¶: Sets a trait notification dispatch handler.

setup_source(fname)¶: Given a VTK file name fname, this creates a mayavi2 reader for it and adds it to the pipeline. It returns the reader created.

trait_added¶: An event fired when a new trait is dynamically added to the object. The value is the name of the trait that was added.

trait_items_event(name, event_object, event_trait)¶

Fire an items event for changes to a Traits collection.

Parameters:

name (str) – Name of the item trait for which an event is being fired. (The name will usually end in ‘_items’.)
event_object (object) – Object of type TraitListEvent, TraitDictEvent or TraitSetEvent describing the changes to the underlying collection trait value.
event_trait (CTrait) – The items trait, of trait type Event.

trait_modified¶: An event that can be fired to indicate that the state of the object has been modified.

trait_property_changed(name, old_value[, new_value])¶

Call notifiers when a trait property value is explicitly changed.

Calls trait and object notifiers for a property value change.

Parameters:

name (str) – Name of the trait whose value has changed
old_value (any) – Old value for this trait.
new_value (any, optional) – New value for this trait. If the new value is not provided, it’s looked up on the object.

classmethod trait_subclasses(all=False)¶

Returns a list of the immediate (or all) subclasses of this class.

Parameters:: all (bool) – Indicates whether to return all subclasses of this class. If False, only immediate subclasses are returned.

traits_init()¶

Perform any final object initialization needed.

For the CHasTraits base class, this method currently does nothing.

traits_inited()¶

Get the initialization state of this object.

Returns:: initialized – True if the object is initialized, else False.
Return type:: bool

update_pipeline(source)¶: Override this to do something else if needed.

wrappers = {'extended': <class 'traits.trait_notifiers.ExtendedTraitChangeNotifyWrapper'>, 'fast_ui': <class 'traits.trait_notifiers.FastUITraitChangeNotifyWrapper'>, 'new': <class 'traits.trait_notifiers.NewTraitChangeNotifyWrapper'>, 'same': <class 'traits.trait_notifiers.TraitChangeNotifyWrapper'>, 'ui': <class 'traits.trait_notifiers.FastUITraitChangeNotifyWrapper'>}¶: Mapping from dispatch type to notification wrapper class type

examples.cahnHilliard.tanh1D module¶

This example solves the Cahn-Hilliard equation given by,

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

where the free energy functional is given by,

$f = \frac{a^2}{2} \phi^2 (1 - \phi)^2$

The Cahn-Hilliard equation can be rewritten in the following form,

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

The above form of the equation makes the non-linearity part of the diffusion coefficient for the first term on the RHS. This is the correct way to express the equation to FiPy.

We solve the problem on a 1D mesh

>>> from fipy import CellVariable, Grid1D, NthOrderBoundaryCondition, DiffusionTerm, TransientTerm, LinearLUSolver, DefaultSolver, Viewer
>>> from fipy.tools import numerix

>>> L = 40.
>>> nx = 1000
>>> dx = L / nx
>>> mesh = Grid1D(dx=dx, nx=nx)

and create the solution variable

>>> var = CellVariable(
...     name="phase field",
...     mesh=mesh,
...     value=1.)

The boundary conditions for this problem are

$\left. \begin{aligned} \phi &= \frac{1}{2} \\ \frac{\partial^2 \phi}{\partial x^2} &= 0 \end{aligned} \right\} \qquad \text{on $x = 0$}$

and

$\left. \begin{aligned} \phi &= 1 \\ \frac{\partial^2 \phi}{\partial x^2} &= 0 \end{aligned} \right\} \qquad \text{on $x = L$}$

or

>>> BCs = (
...     NthOrderBoundaryCondition(faces=mesh.facesLeft, value=0, order=2),
...     NthOrderBoundaryCondition(faces=mesh.facesRight, value=0, order=2))

>>> var.constrain(1, mesh.facesRight)
>>> var.constrain(.5, mesh.facesLeft)

Using

>>> asq = 1.0
>>> epsilon = 1
>>> diffusionCoeff = 1

we create the Cahn-Hilliard equation:

>>> faceVar = var.arithmeticFaceValue
>>> freeEnergyDoubleDerivative = asq * ( 1 - 6 * faceVar * (1 - faceVar))

>>> diffTerm2 = DiffusionTerm(
...     coeff=diffusionCoeff * freeEnergyDoubleDerivative)
>>> diffTerm4 = DiffusionTerm(coeff=(diffusionCoeff, epsilon**2))
>>> eqch = TransientTerm() == diffTerm2 - diffTerm4

>>> import fipy.solvers.solver
>>> if fipy.solvers.solver in ['pysparse', 'pyamgx']:
...     solver = LinearLUSolver(tolerance=1e-15, iterations=100)
... else:
...     solver = DefaultSolver()

The solution to this 1D problem over an infinite domain is given by,

$\phi(x) = \frac{1}{1 + \exp{\left(-\frac{a}{\epsilon} x \right)}}$

or

>>> a = numerix.sqrt(asq)
>>> answer = 1 / (1 + numerix.exp(-a * (mesh.cellCenters[0]) / epsilon))

If we are running interactively, we create a viewer to see the results

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=var, datamin=0., datamax=1.0)
...     viewer.plot()

We iterate the solution to equilibrium and, if we are running interactively, we update the display and output data about the progression of the solution

>>> dexp=-5
>>> from builtins import range
>>> for step in range(100):
...     dt = numerix.exp(dexp)
...     dt = min(10, dt)
...     dexp += 0.5
...     eqch.solve(var=var, boundaryConditions=BCs, solver=solver, dt=dt)
...     if __name__ == '__main__':
...         diff = abs(answer - numerix.array(var))
...         maxarg = numerix.argmax(diff)
...         print('maximum error: {}'.format(diff[maxarg]))
...         print('element id:', maxarg)
...         print('value at element {} is {}'.format(maxarg, var[maxarg]))
...         print('solution value: {}'.format(answer[maxarg]))
...
...         viewer.plot()

We compare the analytical solution with the numerical result,

>>> print(var.allclose(answer, atol=1e-4))
1

examples.cahnHilliard.test module¶

Last updated on Jun 27, 2023. Created using Sphinx 6.2.1.