gistfile1.py

import fipy as fp

Lx1 = 2e-2
Lx2 = 2e-2
Lx3 = 2e-2
Ly = 3e-3
Lx = Lx1 + Lx2 + Lx3
nx = 20
ny = 10
dx = Lx/nx
dy = Ly/ny
mesh = fp.Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
phi = fp.CellVariable(name="elec. potential", mesh=mesh, value=0.)
c_m = fp.CellVariable(name="conc.", mesh=mesh, value=1.0)

V = 0.0007 # V
c0 = 1. # part./m^3
c_m.setValue(c0)
D_m = 2.03e-9 # m^2/s
D_p = 1.33e-9 # m^2/s
Xfc, Yfc = mesh.faceCenters

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)

print "phi con bot"
phi.faceGrad.constrain((D_m - D_p)/(D_p + D_m) * c_m.faceGrad / c_m.faceValue,
        where=mesh.facesBottom)
print "phi con top"
phi.faceGrad.constrain((D_m - D_p)/(D_p + D_m) * c_m.faceGrad / c_m.faceValue,
#phi.faceGrad.constrain(2.0 * c_m.faceGrad,
        where=mesh.facesTop)
print "cm"
c_m.faceGrad.constrain( c_m.faceValue * phi.faceGrad, mesh.facesBottom)
c_m.faceGrad.constrain( c_m.faceValue * phi.faceGrad, mesh.facesTop)

eq = eq1 & eq2

res = 1e10
nIt = 0
while res > 1e-3:
    res = eq.sweep()
    print "Iteration:", nIt, "Residual:", res
    nIt += 1
print "{num} iterations".format(num=nIt)