Internal zero-flux boundary condition in 2D advection-diffusion equation

[taoyun]

I wish to solve the transient solutions of a Fokker-Planck equation on a 2-dimensional spatial domain where there exists an internal rectangular "barrier" of which the boundaries are zero-flux. I've referred to the official documentation (https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#internal-fixed-value) as well as numerous forum posts, all without success. I've set up the following simple example (using FiPy 3.1):

import numpy as np
import fipy
from fipy import numeric
import matplotlib.pyplot as plt
from matplotlib.mlab import bivariate_normal
from matplotlib.patches import Rectangle

fig =plt.figure()
ax = plt.subplot()
plt.show()

# generate mesh
l = 6.
nx = 301 
ny = nx
dx = l/nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]] 

convection = VanLeerConvectionTerm

# reshaping solution for visualization
def phi_reshape(args):
    p1 = args.value.reshape(mesh.shape)
    return p1

# initial condition
init = bivariate_normal(mesh.x, mesh.y, 0.1, 0.1, -2., 2.)
phi = CellVariable(name="solution variable", mesh=mesh, value=init)

# exterior boundary condition
BCs = (FixedFlux(faces=mesh.getFacesRight(),  value=0.),
       FixedFlux(faces=mesh.getFacesLeft(),   value=0.),
       FixedFlux(faces=mesh.getFacesTop(),    value=0.),
       FixedFlux(faces=mesh.getFacesBottom(), value=0.))

# base parameters
d0 = 0.1
c0 = 0.8

# internal rectangle
X, Y = mesh.faceCenters()
mask = ((X > -1.7) & (X < -1.5) & (Y > 1.7))
d = FaceVariable(mesh=mesh, value=d0)
d[mask] = 0.
c = FaceVariable(mesh=mesh, value=c0)
c[mask] = 0.

largeValue = 1e+10
gradient = 0.

eqX = (TransientTerm(var=phi) == DiffusionTerm(coeff=d, var=phi))
+ Diffusion(coeff=largeValue * mask, var=phi)
+ convection(coeff=c, var=phi)
+ convection(coeff=largeValue * mask, var=phi)
- ImplicitSourceTerm(mask * largeValue * gradient * mesh.faceNormals).divergence, var=phi)

# initial plot
p1 = phi_reshape(phi)
img = ax.imshow(p1, origin='lower', extent=[-l/2, l/2, -l/2, l/2])
cb = fig.colorbar(img)
rec = ax.add_patch(Rectangle((-1.7, 1.7), 0.2, 2, color='yellow', alpha=0.2))
ax.set_xlim(-2.5, -1.)
ax.set_ylim(1., 2.5)
plt.draw()

# solve
dt = 0.002
steps = 100

for step in range(steps):
    print 'step', step
    eqX.solve(var=phi, boundaryConditions=BCs, dt=dt)
    p1 = phi_reshape(phi)
    img = ax.imshow(p1, origin='lower', extent=[-l/2, l/2, -l/2, l/2])
    cb.update_bruteforce(img)
    
    ax.set_xlim(-2.5, -1.)
    ax.set_ylim(1., 2.5)
    plt.draw()

I cannot resolve a number of major issues with eqX.

1. The use of two DiffusionTerm alone (i.e., without the convention and the implicit source terms) appears to generate broadcasting errors.
2. I could not compute the divergence in the implicit source term due to its coefficient being treated as just a numpy array.
3. There appears to be no way to incorporate fixed gradient in the convection term.

Any help to enable this model would be helpful to a presumably large number of FiPy users.