This script illustrates a failure with the finite volume method for an elliptical problem on a non-uniform mesh. For more information see the question and answer here, http://scicomp.stackexchange.com/questions/8577/peculiar-error-when-solving-the-poisson-equation-on-a-non-uniform-mesh-1d-only

fipy_nonuniform_poisson_electrostatics.py

from fipy import *
import numpy as np

"""
This script illustrates a failure with the finite volume method for an elliptical problem on a non-uniform mesh.
For more information see the question and answer here,
http://scicomp.stackexchange.com/questions/8577/peculiar-error-when-solving-the-poisson-equation-on-a-non-uniform-mesh-1d-only
"""

def fipy_step_charge():
    nx = 20
    dx  = 0.5
    
    # First generate a uniform mesh.
    # Second, remove one cell and expand it's neigbour 
    # to keep the domain the original size.
    # Removing cell 6 caues no problems, however,
    # removing cell 7 causes the method to fail.
    remove_cell_inx = 6
    width = np.ones(nx-1) * 0.5 
    width[remove_cell_inx] = dx * 2.0
    L = np.sum(width)
    dep = L/10.
    dx = width
    
    mesh = Grid1D(nx=nx-1, dx=dx )
    x = mesh.cellCenters[0]

    potential = CellVariable(mesh=mesh, name='potential', value=0.)
    permittivity = 1

    electrons = CellVariable(mesh=mesh, name='e-')
    electrons.setValue(0)
    electrons.setValue(1., where=x > (L / 2. + dep))
    electrons.valence = -1

    holes = CellVariable(mesh=mesh, name='h+')
    holes.setValue(0)
    holes.setValue(1., where=x < (L / 2. - dep))
    holes.valence = 1

    idonors = CellVariable(mesh=mesh, name='n*')
    idonors.setValue(0)
    idonors.setValue(1., where=x > (L / 2.))
    idonors.valence = 1

    iacceptors = CellVariable(mesh=mesh, name='p*')
    iacceptors.setValue(0)
    iacceptors.setValue(1., where=x < (L / 2.))
    iacceptors.valence = -1

    charge = electrons * electrons.valence + \
             holes * holes.valence +         \
             idonors * idonors.valence +     \
             iacceptors * iacceptors.valence
             
    charge.name = "charge"
    potential.equation = (DiffusionTerm(coeff = permittivity) + charge == 0 )
    potential.constrain(0., mesh.facesLeft)
    potential.faceGrad.constrain(0., mesh.facesRight)
    potential.equation.solve(var=potential)
    
    return (x, np.array(potential), np.array(charge))


if __name__ == '__main__':
    x, phi, charge = fipy_step_charge()
    import pylab
    pylab.plot(x, phi, "o-", label="potential")
    pylab.plot(x, charge, "r--", label="input charge")
    pylab.legend()
    pylab.show()
    print x, phi