The free energy is given by,
$$ f_0\left[ c \left( \vec{r} \right) \right] = - \frac{A}{2} \left(c - c_m\right)^2 + \frac{B}{4} \left(c - c_m\right)^4 + \frac{c_{\alpha}}{4} \left(c - c_{\alpha} \right)^4 + \frac{c_{\beta}}{4} \left(c - c_{\beta} \right)^4 $$In FiPy we write the evolution equation as
$$ \frac{\partial c}{\partial t} = \nabla \cdot \left[ D \left( c \right) \left( \frac{ \partial^2 f_0 }{ \partial c^2} \nabla c - \kappa \nabla \nabla^2 c \right) \right] $$Let's start by calculating $ \frac{ \partial^2 f_0 }{ \partial c^2} $ using sympy. It's easy for this case, but useful in the general case for taking care of difficult book keeping in phase field problems.
In [1]:
%matplotlib inline
import sympy
import fipy as fp
import numpy as np
In [2]:
A, c, c_m, B, c_alpha, c_beta = sympy.symbols("A c_var c_m B c_alpha c_beta")
In [3]:
f_0 = - A / 2 * (c - c_m)**2 + B / 4 * (c - c_m)**4 + c_alpha / 4 * (c - c_alpha)**4 + c_beta / 4 * (c - c_beta)**4
In [4]:
print f_0
In [5]:
sympy.diff(f_0, c, 2)
Out[5]:
The first step in implementing any problem in FiPy is to define the mesh. For Problem 1a the solution domain is just a square domain, but the boundary conditions are periodic, so a PeriodicGrid2D object is used. No other boundary conditions are required.
In [6]:
mesh = fp.Grid2D(nx=50, ny=50, dx=0.5, dy=0.5)
The next step is to define the parameters and create a solution variable.
In [7]:
c_alpha = 0.05
c_beta = 0.95
A = 2.0
kappa = 2.0
c_m = (c_alpha + c_beta) / 2.
B = A / (c_alpha - c_m)**2
D = D_alpha = D_beta = 2. / (c_beta - c_alpha)
c_0 = 0.45
q = np.sqrt((2., 3.))
epsilon = 0.01
c_var = fp.CellVariable(mesh=mesh, name=r"$c$", hasOld=True)
Now we need to define the initial conditions given by,
Set $c\left(\vec{r}, t\right)$ such that
$$ c\left(\vec{r}, 0\right) = \bar{c}_0 + \epsilon \cos \left( \vec{q} \cdot \vec{r} \right) $$
In [8]:
r = np.array((mesh.x, mesh.y))
c_var[:] = c_0 + epsilon * np.cos((q[:, None] * r).sum(0))
viewer = fp.Viewer(c_var)
To define the equation with FiPy first define f_0 in terms of FiPy. Recall f_0 from above calculated using Sympy. Here we use the string representation and set it equal to f_0_var using the exec command.
In [9]:
out = sympy.diff(f_0, c, 2)
In [10]:
exec "f_0_var = " + repr(out)
In [11]:
#f_0_var = -A + 3*B*(c_var - c_m)**2 + 3*c_alpha*(c_var - c_alpha)**2 + 3*c_beta*(c_var - c_beta)**2
f_0_var
Out[11]:
In [12]:
eqn = fp.TransientTerm(coeff=1.) == fp.DiffusionTerm(D * f_0_var) - fp.DiffusionTerm((D, kappa))
eqn
Out[12]:
To solve the equation a simple time stepping scheme is used which is decreased or increased based on whether the residual decreases or increases. A time step is recalculated if the required tolerance is not reached.
In [13]:
elapsed = 0.0
steps = 0
dt = 0.01
total_sweeps = 2
tolerance = 1e-1
total_steps = 100
In [14]:
c_var[:] = c_0 + epsilon * np.cos((q[:, None] * r).sum(0))
c_var.updateOld()
from fipy.solvers.pysparse import LinearLUSolver as Solver
solver = Solver()
while steps < total_steps:
res0 = eqn.sweep(c_var, dt=dt, solver=solver)
for sweeps in range(total_sweeps):
res = eqn.sweep(c_var, dt=dt, solver=solver)
if res < res0 * tolerance:
steps += 1
elapsed += dt
dt *= 1.1
c_var.updateOld()
else:
dt *= 0.8
c_var[:] = c_var.old
viewer.plot()
print 'elapsed_time:',elapsed
The following cell will dumpy a file called fipy_hackathon1b.py to the local file system to be run. The images are saved out at each time step.
In [17]:
%%writefile fipy_hackathon_1b.py
import fipy as fp
import numpy as np
mesh = fp.Grid2D(nx=400, ny=400, dx=0.5, dy=0.5)
c_alpha = 0.05
c_beta = 0.95
A = 2.0
kappa = 2.0
c_m = (c_alpha + c_beta) / 2.
B = A / (c_alpha - c_m)**2
D = D_alpha = D_beta = 2. / (c_beta - c_alpha)
c_0 = 0.45
q = np.sqrt((2., 3.))
epsilon = 0.01
c_var = fp.CellVariable(mesh=mesh, name=r"$c$", hasOld=True)
r = np.array((mesh.x, mesh.y))
c_var[:] = c_0 + epsilon * np.cos((q[:, None] * r).sum(0))
f_0_var = -A + 3*B*(c_var - c_m)**2 + 3*c_alpha*(c_var - c_alpha)**2 + 3*c_beta*(c_var - c_beta)**2
eqn = fp.TransientTerm(coeff=1.) == fp.DiffusionTerm(D * f_0_var) - fp.DiffusionTerm((D, kappa))
elapsed = 0.0
steps = 0
dt = 0.01
total_sweeps = 2
tolerance = 1e-1
total_steps = 300
c_var[:] = c_0 + epsilon * np.cos((q[:, None] * r).sum(0))
c_var.updateOld()
from fipy.solvers.pysparse import LinearLUSolver as Solver
solver = Solver()
viewer = fp.Viewer(c_var)
while steps < total_steps:
res0 = eqn.sweep(c_var, dt=dt, solver=solver)
for sweeps in range(total_sweeps):
res = eqn.sweep(c_var, dt=dt, solver=solver)
print ' '
print 'steps',steps
print 'res',res
print 'sweeps',sweeps
print 'dt',dt
if res < res0 * tolerance:
steps += 1
elapsed += dt
dt *= 1.1
if steps % 1 == 0:
viewer.plot('image{0}.png'.format(steps))
c_var.updateOld()
else:
dt *= 0.8
c_var[:] = c_var.old
The movie of the evolution for 300 steps.
The movie was generated with the output files of the form image*.png using the following commands,
$ rename 's/\d+/sprintf("%05d",$&)/e' image*
$ ffmpeg -f image2 -r 6 -i 'image%05d.png' output.mp4
In [18]:
from IPython.display import YouTubeVideo
scale = 1.5
YouTubeVideo('S-EYiO0EltM', width=420 * scale, height=315 * scale, rel=0)
Out[18]:
In [ ]: