from fipy import *
L=10.0e-9
nx=400
ipt=nx/100

dx=L/nx
mesh = Grid1D(nx=nx, dx=dx)
phi = CellVariable(name=r"$\phi$", mesh=mesh)

class CalculatedDiffusionCoefficient(FaceVariable):
    def __init__(self, D, PHI):
        FaceVariable.__init__(self, mesh=PHI.mesh)
        self.PHI = self._requires(PHI)
        self.D = self._requires(D)

    def _calcValue(self):
        arr = numerix.empty(self.PHI.shape)
        for i in xrange(self.PHI.mesh.numberOfFaces):
            arr[i]=self.D*(13077.922*(1/(self.PHI[i]-self.PHI[i]**2))-96000.0)
        print "diffusion coefficient....","        PHI value at",ipt,"=",self.PHI[ipt]
        return arr



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=CalculatedDiffusionCoefficient(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

    print " "
    print "Solving for dt=",dt," PHI value at",ipt,"=",PHI[ipt].value
    eq.solve(phi, dt=dt)
    print "Solved  for dt=",dt," PHI value at",ipt,"=",PHI[ipt].value
    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!')