FiPy electrolyte, transient

gistfile1.py

import fipy as fp

Lx = 1.
nx = 1500
dx = Lx/nx
mesh = fp.Grid1D(nx=nx, dx=dx)
phi = fp.CellVariable(name="elec. potential", mesh=mesh, value=0.,
        hasOld=True)
c_m = fp.CellVariable(name="conc.", mesh=mesh, value=1.0, hasOld=True)

V = -2.5e-0
c0 = 1.
c_m.setValue(c0)
D_m = 2.03e-9 # m^2/s
D_p = 1.33e-9 # m^2/s
D_p = D_p/D_m
D_m = 1.

eq1 = (0 == fp.TransientTerm(coeff=1., var=c_m)
        + fp.DiffusionTerm(coeff=-D_m, var=c_m)
        + fp.DiffusionTerm(coeff=D_m*c_m.faceValue, var=phi))
eq2 = (0 == 0.
        + fp.DiffusionTerm(coeff=-(D_p - D_m), var=c_m)
        + fp.DiffusionTerm(coeff=-(D_p + D_m)*c_m.faceValue, var=phi))

phi.constrain(0., mesh.facesLeft)
phi.constrain(V, mesh.facesRight)
c_m.constrain(c0, mesh.facesLeft)
c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight)

Xcc = mesh.cellCenters[0]
J = 1 - fp.numerix.exp(V)
c_m_analyt = fp.CellVariable(name="analyt. conc.", mesh=mesh)
c_m_analyt.setValue(1 - J*Xcc)
phi_analyt = fp.CellVariable(name="analyt. phi", mesh=mesh)
phi_analyt.setValue(fp.numerix.log(1-J*Xcc))

eq = eq1 & eq2

viewer = fp.Viewer(vars=(c_m, phi))
dt = 5e-4
tf = 8.
nt = tf/dt
t = 0.
success = True
nItMax = 200
while t < tf:
    if success:
        c_m.updateOld()
        phi.updateOld()
    res = 1e10
    nIt = 0
    while res > 5e-2 and nIt < nItMax:
        res = eq.sweep(dt=dt,
                solver=fp.LinearGMRESSolver(),
                )
        # Print every 10th sweep
        if not nIt % 10:
            print "Iteration:", nIt, "Residual:", res
        nIt += 1
    if nIt < 5:
        dt *= 1.01
        success == True
    elif nIt > nItMax:
        dt *= 0.5
        success = False
    t += dt
    print "\n time: {t}, dt: {dt}\n".format(t=t, dt=dt)
    viewer.plot()

anstol = 1e-3
print c_m.allclose(c_m_analyt, rtol=anstol, atol=anstol)
print phi.allclose(phi_analyt, rtol=anstol, atol=anstol)