FiPy quasi-neutral 1D electrolyte

gistfile1.py

import fipy as fp

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

#V = -6e-1
V = -4e-1
#V = 4.0e0
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.DiffusionTerm(coeff=-D_m, var=c_m) +
            fp.DiffusionTerm(coeff=D_m*c_m, var=phi))
eq2 = (0 == fp.DiffusionTerm(coeff=-(D_p - D_m), var=c_m) +
            fp.DiffusionTerm(coeff=-(D_p + D_m)*c_m, 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

res = 1e10
nIt = 0
while res > 1e-6:
    res = eq.sweep()
    print "Iteration:", nIt, "Residual:", res
    nIt += 1

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