def function(D,PHI):
    from copy import deepcopy
    lsize=numerix.size(PHI)
    arr=deepcopy(PHI)
    arr=D*(13077.922*(1/(PHI-PHI**2))-96000.0)
    return arr



from fipy import *
L=10.0e-9
nx=400
dx=L/nx
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(name=r"$\phi$", mesh=mesh)

noice = GaussianNoiseVariable(mesh=mesh,mean=0.5,variance=0.001)
pi=3.14159265359
for i in range(nx):
    z=i*dx/L*2*pi*30
    y=0.5+0.1*numerix.sin(z)
    noice[i] = y

phi.setValue(noice)

PHI = phi.getArithmeticFaceValue()
D = 1e-23

viewer = Viewer(vars=(phi,), datamin=0., datamax=1.)
dexp = -5
elapsed = 0.
dt=1.0e-8
x_new=[]
x_old=[]

eq = TransientTerm() == DiffusionTerm(coeff=function(D,PHI))- DiffusionTerm(coeff=(D, 7.8e-16))

while elapsed < 1.2:
    dt = min(100,numerix.exp(dexp))
    x_old=PHI.all()

    elapsed += dt
    dexp += 0.01

    eq.solve(phi, dt=dt)
    x_new = PHI.all()
    if ((x_new-x_old)/x_old)>0.01:
        dt=dt/2

    viewer.plot()
    print 'dt: ',dt,elapsed


raw_input('Done!')