Inaccuracy when applying internal fixed gradient: evident leakage through zero-flux boundaries

Question

I implemented zero-flux internal boundary condition in 2d space using the official method, i.e., through the use of an implicit source term (https://www.ctcms.nist.gov/fipy/documentation/USAGE.html#internal-fixed-value). In the simple diffusion example below, there exists an internal rectangular "barrier" of which the boundaries are reflective. However, it is clear from the resultant plots (short animation attached) that the solution "leaks" across these supposedly shielded internal boundaries. Is this a limitation of FiPy or might there be a workaround to this issue?

import numpy as np
from fipy import *
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt

plt.rc('text', usetex=True)

fig =plt.figure()
txt = fig.suptitle('place holder', fontsize=14)
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]]

# initial condition
x_init, y_init = -2., 2.
x_var, y_var = 0.1, 0.1
distr = bivariate_normal(mesh.x, mesh.y, x_var, y_var, x_init, y_init)
phi = CellVariable(name="solution variable", mesh=mesh, value=distr)

# function to reshape solution for plotting
def data_process(arg):
    p1 = arg.value.reshape(mesh.shape)
    return p1

# baseline diffusion coefficient
d0 = 0.1

# positional variables
xy_cell = mesh.cellCenters()
xy_face = FaceVariable(mesh=mesh, value=mesh.faceCenters(), rank=1)

# internal zero-flux boundaries
xx, yy = xy_face
mask = ((xx > -1.7) & (xx < -1.5) & (yy > 1.7))
mask = FaceVariable(mesh, value=mask)
d = FaceVariable(mesh=mesh, value=d0)
d[mask.value] = 0.

largeValue = 1e+10

# main equation
eqX = (TransientTerm(var=phi) == DiffusionTerm(coeff=d, var=phi) + \
                                 DiffusionTerm(coeff=largeValue * mask, var=phi) - \
                                 ImplicitSourceTerm(coeff=(mask * largeValue * 0. * mesh.faceNormals).divergence, var=phi))

""" SIMULATION """
# time steps
dt = 0.004
steps = 100

# plot at initial time
p1 = data_process(phi)
txt.set_text(r'Time Step 0')
img = ax.imshow(p1, origin='lower', extent=[-l/2,l/2,-l/2,l/2])
cb = fig.colorbar(img)
plt.draw()

# solving transient solutions
for step in range(steps):
    print 'step', step
    eqX.solve(var=phi, dt=dt)
    
    p1 = data_process(phi)
    txt.set_text(r'Time Step %s'%(step+1))
    
    ax.cla()
    img = ax.imshow(p1, origin='lower', extent=[-l/2,l/2,-l/2,l/2])
    plt.draw()