Space dependent diffusion in fipy

Question

I'm using fipy to model the linearized Poisson-Boltzmann equation, which is essentially

I'll assume I can model f(x) as a boundary condition. If epsilon(x) is a constant, fipy can handle this:

phi = CellVariable(mesh)

dielectric_solvent = 80.0
dielectric_inner   = 4.0

LHS = (DiffusionTerm(coeff = dielectric_solvent))
RHS = phi
eq = LHS == RHS

dr = np.linalg.norm(mesh.faceCenters, axis=0)
mask = (dr<.5) * mesh.exteriorFaces
phi.constrain(1, mask)

mask = (dr>.5) * mesh.exteriorFaces
phi.constrain(0, mask)   

sol = eq.solve(var=phi)

giving:

The complete minimal example is posted as a gist, to keep things short this is the relevant portion.

What I'd like to do is let epsilon(x) vary as a function in space, but DiffusionTerm can only take a constant. How can I implement the spatially varying dielectric term?

Answer 1

Any coefficient in FiPy can be a function of space. For example, you can set the diffusion coefficient as follows,

diffusion_coefficient = dielectric_solvent * ((mesh.x > -0.5) & (mesh.x < 0.5))

and then use that directly in your equation

LHS = (DiffusionTerm(coeff = diffusion_coefficient))

Just define your spatially varying functions using the cell centers, mesh.x and mesh.y.

Another pointer, it might be better to change your equation to

eq = TransientTerm() == DiffusionTerm(diffusion_coefficient) - ImplicitSourceTerm(phi)

so that

• the variable phi is solved implicitly in one sweep

• adding in a TransientTerm stabilizes the problem for me, just use a very large time step to approximate a steady state problem

See my comment on your gist for the changes.