How do you specifiy a Neumann (fixed flux normal to face) boundary condition in Fipy?
<p>How do I explicitly set the flux normal to a boundary face in a fipy mesh to be a specific value, without constraining components of flux within the face?</p>

<p>A Neumann boundary condition can be specified as: (1) fixed component of flux normal to a boundary face, or (2) as a complete specification of flux at the face. The default fipy condition is the former (value = 0), but the explicit method (faceGrad.constrain) is the latter. The issue can be understood by adding the following code to the end of the fipy <a href="http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh20x20.html" rel="nofollow">diffusion.mesh20x20</a> example and noting different face gradient results. </p>

<pre><code>facesNeumann = mesh.exteriorFaces &amp; ~facesTopLeft &amp; ~facesBottomRight
print 'grad(phi) with default Neumann BC...'
print phi.faceGrad.value.T[facesNeumann.value]
phi.faceGrad.constrain(0.,facesNeumann)
DiffusionTerm().solve(var=phi)
print 'and with explicit Neumann BC...'
print phi.faceGrad.value.T[facesNeumann.value]
</code></pre>

<p>The discussion on fixed flux boundary conditions does answer my question, but I had not understood the documentation, which is very brief. Below is a worked example that illustrates how to apply this in a simple 2D case, similar to the mesh20x20 example.</p>

<pre><code>import matplotlib.pyplot as plt
from fipy import *

nx = 20
ny = nx
dx = 1.
dy = dx
L = dx * nx
msh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

xFace, yFace = msh.faceCenters
xCell,yCell = msh.cellCenters    
phi = CellVariable(name = "solution variable",
                   mesh = msh,
                   value = 0.)    
D = FaceVariable(name='diffusion coefficient',mesh=msh,value=1.)

# Dirichlet condition on top face
valueTop = FaceVariable(name='valueTop',mesh=msh,value=xFace*0.1-1)
phi.constrain(valueTop, msh.facesTop)

# Fixed flux (Neumann) on base
D.constrain(0., msh.facesBottom)
fluxBottom = -0.05

eq = DiffusionTerm(coeff=D) + (msh.facesBottom * fluxBottom).divergence
eq.solve(var=phi)

# confirm only y component of gradient is zero
print phi.faceGrad.value.T[msh.facesBottom.value]

plt.scatter(phi.value, yCell)
plt.show()

viewer = Viewer(vars=phi, datamin=-1., datamax=1.)
viewer.plot()
</code></pre>

python mesh pde fipy