CODEWITHSUNDEEP
pythonmeshpdefipy
Rupert
Rupert

Reputation: 31

How do you specifiy a Neumann (fixed flux normal to face) boundary condition in Fipy?

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?

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 diffusion.mesh20x20 example and noting different face gradient results. 

facesNeumann = mesh.exteriorFaces & ~facesTopLeft & ~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]

Upvotes: 2

Views: 1794

Answers (3)

Rupert
Rupert

Reputation: 31

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.

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()

Upvotes: 1

wd15
wd15

Reputation: 1068

In terms of just specifying the flux normal to a boundary for a variable independent of solving an equation, there does not appear to be any way to do that. The syntax might be phi.faceGrad[0].constrain(...) for example, but that currently doesn't work in FiPy. It might also be difficult to implement for arbitrary oriented faces.

However, for practical purposes, the component tangential to a boundary is not used when solving an equation in FiPy, only the normal component has any influence on the solution. For example, take the following code,

import fipy as fp
mesh = fp.Grid2D(nx=2, ny=1)
var = fp.CellVariable(mesh=mesh)
var.constrain(1, mesh.facesLeft)
var.constrain(0, mesh.facesRight)
#var.faceGrad.constrain(0, mesh.facesTop)
fp.DiffusionTerm().solve(var)
print 'face gradient on top plane:',var.faceGrad[0, mesh.facesTop.value]
print 'variable value:',var

This gives an output of

face gradient on top plane: [-0.5 -0.5]
variable value: [ 0.75  0.25]

The answer is correct, but the top face gradient is -0.5. However, when the line #var.faceGrad.constrain(0, mesh.facesTop) is uncommented, the result is

face gradient on top plane: [ 0.  0.]
variable value: [ 0.75  0.25]

The tangential face gradient is now 0, but the answer is the same. The point is that setting the face gradient tangentially (via .constrain) has no influence on the solution.

Upvotes: 2

jeguyer
jeguyer

Reputation: 2484

Please see the discussion on fixed flux boundary conditions. Basically, you add a source containing the the divergence of the desired boundary flux to FiPy's default no-flux condition.

Upvotes: 2

Related Questions