#!/usr/bin/env python

#

nx = 100
ny = 100
dx = 1e-6 # cm
dy  = 1e-6
L = nx * dx
from fipy.meshes.grid2D import Grid2D
mesh = Grid2D(dx=dx, nx=nx, dy=dy , ny=ny)
#

#
from fipy.tools.numerix import random
from fipy.variables.cellVariable import CellVariable

C_li = CellVariable(name="Lithium Ion concentration", mesh=mesh, value=1.0, hasOld=1)
POT = CellVariable(name="Potential", mesh=mesh, value=1.0, hasOld=1)
# ..




from fipy.tools.numerix import log

epsilon= (8.85419e-14)*(8.0) #f/cm   epsilon_not*8
epsilon_prime=0.0
#
Vm = 12.9981 # cm**3 / mol








#




Mphi = 10000.0
D_electrons=0.00001 # cm**2/s
D_li=0.00001 # cm**2/s
#
M_electrons=0.00001  #??
M_li=0.00001  #??

#
from fipy.terms.transientTerm import TransientTerm
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.implicitSourceTerm import ImplicitSourceTerm
from fipy.terms.diffusionTerm import DiffusionTerm
from fipy.terms.powerLawConvectionTerm import PowerLawConvectionTerm
C_li = CellVariable(name="concentration_lithium_ions", mesh=mesh, hasOld=1)

Rho= 96485.3 * C_li

amp=L/16.0
nwave=3.0
pi=3.1415926535898
from fipy.tools.numerix import sin
C_li.setValue(1.0)
C_li.setValue(0.25, where=mesh.getCellCenters()[...,0] > L/7.0+amp*sin(nwave*2.0*pi*mesh.getCellCenters()[...,1]/L))



#
from fipy.boundaryConditions.fixedValue import FixedValue
BCS=(FixedValue(mesh.getFacesRight(),0),
     FixedValue(mesh.getFacesLeft(),0))

  
POT_Eq = Rho + ImplicitDiffusionTerm(epsilon)
#

#

#
amp=L/16.0
nwave=3.0
pi=3.1415926535898
from fipy.tools.numerix import sin



#

# ..
#
from fipy.tools.numerix import exp








# ..
#
if __name__ == '__main__':
    from fipy import viewers
  
                         
    POT_viewer =  viewers.make(vars=(POT),  limits={'datamin': 0.0, 'datamax': 1.0})

    POT_viewer.plot()



from fipy.solvers.linearLUSolver import LinearLUSolver
solver = LinearLUSolver(tolerance=1e-10)




 

#
dt = 1.e-6
t=0
#
for i in range(100):

    print i
 
    POT.updateOld()
    t=t+dt

    POTRes = 1e+10
    while  POTRes > 1e-3:
        
        POTRes =  POT_Eq.sweep(var=POT,solver=solver,boundaryConditions = BCS)
    if __name__ == '__main__':

        POT_viewer.plot()
        print t
#
if __name__ == '__main__': 
    raw_input("moving phi field")
#