# Table of Contents
* [1b. Fixed flux spinodal decomposition on a square domain](#1b.-Fixed-flux-spinodal-decomposition-on-a-square-domain)
	* [Use Binder For Live Examples](#Use-Binder-For-Live-Examples)
	* [Define $f_0$](#Define-$f_0$)
	* [Define the Equation](#Define-the-Equation)
	* [Solve the Equation](#Solve-the-Equation)
	* [Run the Example Locally](#Run-the-Example-Locally)
	* [Movie of Evolution](#Movie-of-Evolution)


# 1b. Fixed flux spinodal decomposition on a square domain

## Use Binder For Live Examples

[](http://mybinder.org/repo/wd15/fipy-hackathon1)

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.


```python
%matplotlib inline

import sympy
import fipy as fp
import numpy as np
```


```python
A, c, c_m, B, c_alpha, c_beta = sympy.symbols("A c_var c_m B c_alpha c_beta")
```


```python
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
```


```python
print f_0
```

    -A*(-c_m + c_var)**2/2 + B*(-c_m + c_var)**4/4 + c_alpha*(-c_alpha + c_var)**4/4 + c_beta*(-c_beta + c_var)**4/4



```python
sympy.diff(f_0, c, 2)
```




    -A + 3*B*(c_m - c_var)**2 + 3*c_alpha*(c_alpha - c_var)**2 + 3*c_beta*(c_beta - c_var)**2



The first step in implementing any problem in FiPy is to define the mesh. For [Problem 1a]({{ site.baseurl }}/hackathon1/problems.ipynb/#a.-Square-Periodic) 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.


```python
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.


```python
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) $$


```python
r = np.array((mesh.x, mesh.y))
c_var[:] = c_0 + epsilon * np.cos((q[:, None] * r).sum(0))
viewer = fp.Viewer(c_var)
```

## Define $f_0$

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.


```python
out = sympy.diff(f_0, c, 2)
```


```python
exec "f_0_var = " + repr(out)
```


```python
#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
```




    (((((pow((0.5 - $c$), 2)) * 29.6296296296) + -2.0) + ((pow((0.05 - $c$), 2)) * 0.15)) + ((pow((0.95 - $c$), 2)) * 2.85))



## Define the Equation


```python
eqn = fp.TransientTerm(coeff=1.) == fp.DiffusionTerm(D * f_0_var) - fp.DiffusionTerm((D, kappa))
eqn
```




    (TransientTerm(coeff=1.0) + (DiffusionTerm(coeff=[-(((((((pow((0.5 - $c$), 2)) * 29.6296296296) + -2.0) + ((pow((0.05 - $c$), 2)) * 0.15)) + ((pow((0.95 - $c$), 2)) * 2.85)) * 2.22222222222))]) + DiffusionTerm(coeff=[2.2222222222222223, 2.0])))



## Solve the Equation

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.


```python
elapsed = 0.0
steps = 0
dt = 0.01
total_sweeps = 2
tolerance = 1e-1
total_steps = 100
```


```python
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
```

## Run the Example Locally

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.


```python
%%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
```

    Writing fipy_hackathon_1b.py


## Movie of Evolution

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


```python
from IPython.display import YouTubeVideo
scale = 1.5
YouTubeVideo('S-EYiO0EltM', width=420 * scale, height=315 * scale, rel=0)
```











```python

```

d483c531939e410690415afdb7f2b950b265ed5f
ipynb
Jupyter Notebook
hackathon1/fipy/1b.ipynb
dschwen/chimad-phase-field
a7248ee5af896719d8b604faaa2781bdbd8cbe64
MIT
hackathon1/fipy/1b.ipynb
dschwen/chimad-phase-field
a7248ee5af896719d8b604faaa2781bdbd8cbe64
MIT
hackathon1/fipy/1b.ipynb
dschwen/chimad-phase-field
a7248ee5af896719d8b604faaa2781bdbd8cbe64
MIT
Qwen/Qwen-72B
1. YES
2. YES