Dear FiPy users --
This is a simple example of how and why fipy may fail to solve a
steady advection diffusion problem, and how solving the transient
problem can reduce the error. I also found something that was a
surprise -- the "initial" condition of a steady problem can affect
the solution for some solvers.
At the end are two interesting questions for those who want to
understand what FiPy is actually doing.... I admit to being a bit
lost
The equation I am solving is
\Del\dot (D\del psi + u*psi) =0
Where the diffusion D is 1, and u is a vector (1,0) -- so there is
only a flow of speed -1 in the x direction. I solve this equation
on a 10 by 10 grid. The boundary conditions are no normal gradient
on the y=0 and y=10 boundary:
psi_y =0 at y=0 and y=10
For the x boundary, I impose a value of x=1 on the inflow boundary at
x=10
(this is a little tricky -- the way the equation is written, u is
the negative of velocity).
psi=1 at x=10
and a no-normal-gradient condition at x=0.
psi_x=0 at x=0
since all of the domain and boundary is symmetrical with respect to
y, we can assume psi_y=0 is zero everywhere. This reduces the equation to
psi_xx + psi_x =0
The general solution to this equation is
psi=C1+C2*exp(-x)
Where C1 and C2 are constants. For these boundary conditions, C1=1
and C2=0, so psi=1 everywhere.
Now run the code SquareGrid_HomemadeDelaunay and look at figure(3)
-- this is the plot of psi versus x, and you can see that it does
not match the true solution of psi=1 everywhere! Instead, it
appears to be decaying exponential. The blue line is actually just
(1+exp(-x)). What is going on? It appears that small errors in the
boundary condition at x=0 are allowing C2 to not be exactly 0, and
this error is this exponential mode. The error is the artificial
exiting of a correct mode of the interior equation, albeit one that
should not be excited by these BC's.
Interestingly, the size of this error depends on the value the psi
is initially set to. If the line
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
is changed so psi is initially 1, the error goes away entirely; if
it is set to some other value, you get different errors. I do not
know if this is a bug, or just numerical error in a well programmed
solver.
Now if you run SquareGrid_HomemadeDelaunay_transient which implements
psi_t = \Del\dot (D\del psi + u*psi)
you can see that the error in the numerical solution is advected
out of the x=0 boundary, and the solution approaches the true
solution of psi=1 rapidly.
The interesting question is if the formulation of the boundary
condition at x=0 could be altered to less excite the spurious mode?
Also, why does the "initial" condition have any effect on the
steady equation?
Cheers,
Jamie
from pylab import *
from numpy import *
import fipy
from collections import defaultdict
import copy as Copy
def ScipyDelaunay2mesh(points,simplices,debug=False):
'''ScipyDelaunay2mesh(points,simplices):
Creates a FiPy 2Dmesh object from the output of a
scipy.spatial.Delaunay triangulation,
Also returns a list of border faces for the application of
boundary conditions. The list of border conditions is computed
within this routine; it will handle the case where
triangles/simplices have been removed from the output of the
Delaunay() triangulation.
Parameters
----------
points : array
an (N,2) array of the location of N points that made the
triangulation. If Delaunay is called with
tri=scipy.spatial.Delaunay(), points is tri.points
simplices : array
an (M,3) array of the location of M triangles returned by
Delaunay. If Delaunay is called with
tri=scipy.spatial.Delaunay(), simplices is tri.simplices
debug=False : Boolean
if this is set to true, the mesh will be drawn in blue,
boundaries will be drawn in red, boundary point numbers will be
written at vertices, and numbers for triangles will be drawn in
the triangle. This output is useful for finding where triangles
have been removed incorrectly.
Returns
-------
mesh: a FiPy mesh2D object
'''
#Create the arrays that fipy.mesh.mesh2D needs. vertexCoords is
#of shape (2,N), where N is the number of points, and contains the
#coordinates of the vertices. It is equivalent to points in
#content, but reshaped.
print 'Entering ScipyDelaunay2Mesh; Making vertexCoords'
vertexCoords=points.T
#faceVertexID is of shape (2,M) where M is the number of faces/edges
#(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
#terminology). faceVertexID contains the points (as indices into
#vertexCoords) which make up each face. This data can be extracted
#from simplices from Delaunay, which contains the faces of each
#triangle. HOWEVER, simplices contains many repeated faces, since each
#triangle will border on others. Only one copy of each face should be
#inserted into faceVertexID.
print ' Making faceVertexIDs'
faceIn={} #this is a dictionary of all faces that have already been
#inserted into the faceVertexID, with the key a tuple of
#indices, sorted
faceVertexIDs_list=[] #structure to build
for ntri in xrange(simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=simplices[ntri,[0,1]]
elif nface==1:
face=simplices[ntri,[1,2]]
else:
face=simplices[ntri,[2,0]]
faceKey=tuple(sort(face))
if not (faceKey in faceIn):
faceIn[faceKey]=True
faceVertexIDs_list.append(face)
faceVertexIDs=array(faceVertexIDs_list).T
#now create cellFaceIDs of shape (3,P) where P is the number of cells
#in the domain. It contains the faces (as indices into faceVertexIDs)
#that make up each cell.
#first create dictionary with a key made up of a sorted tuple of vertex
#indices which maps from a face to its location in faceVertexIDs
print ' Making cellFaceIDs'
faceMap={}
for nface in xrange(faceVertexIDs.shape[1]):
faceKey=tuple(sort(faceVertexIDs[:,nface]))
faceMap[faceKey]=nface
#now loop over simplices, find each edge/face, and map from points to
#faces
cellFaceIDs=simplices.T*0
for ntri in xrange(simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=simplices[ntri,[0,1]]
elif nface==1:
face=simplices[ntri,[1,2]]
else:
face=simplices[ntri,[2,0]]
faceKey=tuple(sort(face))
cellFaceIDs[nface,ntri]=faceMap[faceKey]
#ok, now instantiate a mesh2D object with this data
print ' Making mesh'
mesh=fipy.meshes.mesh2D.Mesh2D(vertexCoords,faceVertexIDs,cellFaceIDs)
return mesh
def Matrix2Delaunay(xmat,ymat,mask,addOrthog=False):
'''Matrix2Delaunay(xmat,ymat,mask,addOrthog=False):
Creates the elements of a Delaunay triangulation, simplices and
points, from two matrices, containing the x and y positions of the
points, and a mask array which. Only include triangles connecting
points where mask is true.
Parameters
----------
xmat: x location of points, 2D matrix
ymat: y location of points, 2D matrix
mask: boolean array, 2D matrix
mask.shape must be the same as xmat.shape and ymat.shape
mask MUST be so that triangles are able to be formed!
addOrthog=False: if true add points to result which improve
orthogonality of resulting triangles.
Returns
-------
points, the location of vertices as returned by Delaunay
simplices, the points which make up triangles, as returned by Delaunay
pntMat: a matrix of xmat.shape which contains the indices of the
vertices that correspond to points in the matrix. Note -- does
not return points that where added to improve orthogonality,
since they are not in xmat,ymat.
'''
assert xmat.shape==mask.shape, 'xmat, ymat and mask must be same shape'
assert ymat.shape==mask.shape, 'xmat, ymat and mask must be same shape'
assert mask.dtype==type(True), 'mask must be a boolean array'
#make lists of x and y points that are to be included in mesh, and
#pntMat, a matrix that maps from a location in the patrix to the
#index of the points.
xvec=xmat[mask]
yvec=ymat[mask]
pntMat=xmat.astype(type(1))
pntMat[:,:]=-1 #-1 means point not included in mesh
pntMat[mask]=arange(sum(mask))
#Create some of the arrays that Delaunay would create. points is
#of shape (N,2), where N is the number of points, and contains the
#coordinates of the vertices.
print 'Entering Matrix2Delaunay; Making points'
points=zeros([len(xvec),2])
points[:,0]=xvec
points[:,1]=yvec
#now create simplices, an (M,3) array of integers for the M
#triangles that will be created.
print ' making simplices'
if not(addOrthog):
#If the nodes are arranged
#
# 4--3
# | |
# 1--2
#
#create triangles 1->2->4 and 2->3->4
simplices_list=[]
for nx in range(xmat.shape[0]-1):
for ny in range(xmat.shape[1]-1):
#if possible, create 1->2->4
if mask[nx,ny] & mask[nx+1,ny] & mask[nx,ny+1]:
simplices_list.append([pntMat[nx,ny],pntMat[nx+1,ny],pntMat[nx,ny+1]])
#if possible, create 2->3->4
if mask[nx+1,ny] & mask[nx+1,ny+1] & mask[nx,ny+1]:
simplices_list.append([pntMat[nx+1,ny],pntMat[nx+1,ny+1],pntMat[nx,ny+1]])
else:
#If the nodes are arranged
#
# 4--3
# | |
# 1--2
#
#if (1->2) or (2->3) or (3->4) or (4->1) exist, create a new point np
#
# 4----3
# | |
# | np |
# | |
# 1----2
#
#create triangles 1->2->np, 2->3->np, 3->4->np and 4->1->np
#create list from points, same shape as points, to append new points to
points_list=[list(np) for np in list(points)]
simplices_list=[]
for nx in range(xmat.shape[0]-1):
for ny in range(xmat.shape[1]-1):
#only make points if the following sets are not masked
#(1,2,3) or (2,3,4) or (3,4,1) to avoid jaggies on the
#boundary
if ((mask[nx,ny] & mask[nx+1,ny] & mask[nx+1,ny+1]) |
(mask[nx+1,ny] & mask[nx+1,ny+1] & mask[nx,ny+1]) |
(mask[nx+1,ny+1] & mask[nx,ny+1] & mask[nx,ny])):
xnp=0.25*(xmat[nx,ny]+xmat[nx+1,ny]+xmat[nx+1,ny+1]+xmat[nx,ny+1])
ynp=0.25*(ymat[nx,ny]+ymat[nx+1,ny]+ymat[nx+1,ny+1]+ymat[nx,ny+1])
points_list.append([xnp,ynp])
n_newpoint=len(points_list)-1
if (mask[nx,ny] & mask[nx+1,ny]): # add triangle 1->2->np
simplices_list.append([pntMat[nx,ny],pntMat[nx+1,ny],n_newpoint])
if (mask[nx+1,ny] & mask[nx+1,ny+1]): # add triangle 2->3->n
simplices_list.append([pntMat[nx+1,ny],pntMat[nx+1,ny+1],n_newpoint])
if (mask[nx+1,ny+1] & mask[nx,ny+1]): # add triangle 3->4->np
simplices_list.append([pntMat[nx+1,ny+1],pntMat[nx,ny+1],n_newpoint])
if (mask[nx,ny+1] & mask[nx,ny]): # add triangle 4->1->np
simplices_list.append([pntMat[nx,ny+1],pntMat[nx,ny],n_newpoint])
#convert points_list back to an array
points=array(points_list)
#convert list to array
simplices=array(simplices_list)
#what happens if we reverse order of triangles
if False:
for ns in range(simplices.shape[0]):
j0=simplices[ns,0]
j2=simplices[ns,2]
simplices[ns,0]=j2
simplices[ns,2]=j0
return points,simplices,pntMat
if __name__ == "__main__":
#test code by making and drawing simple grid
#make points in square domain of size 1
points=[]
xvec=linspace(0.0,10.0,50)
yvec=linspace(0.0,10.0,50)
xmat,ymat=meshgrid(xvec,yvec)
#make Mask
mask=xmat<Inf
#now make triangulation, and convert to mesh
points,simplices,matPnts=Matrix2Delaunay(xmat,ymat,mask,addOrthog=True)
mesh=ScipyDelaunay2mesh(points,simplices)
#setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
#first, make grid variable
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
figure(1)
clf()
#now plot mesh
triplot(mesh.vertexCoords[0,:],mesh.vertexCoords[1,:],triangles=simplices)
title('simple square mesh')
axis('square')
draw()
show()
# These modules takes the output of scipy.spatial.Delaunay's
# triangulation of data points, and returns both the FiPy mesh for the
# points, and a list of the boundary edges so that boundary conditions
# can be set for the points.
#
# This code handles the case where triangles/simplices have been
# removed from output of the Delaunay triangulation after the
# triangulation. This triangle removal can be useful to create
# interior holes in the domain, or to handle concave boundaries of the
# domain.
#
# The documentation for the routines can be found in there docstrings,
# and examples of there use can be found in __main__
from numpy import *
import fipy
from collections import defaultdict
import copy as Copy
from JMP_Make_Grids import *
if __name__ == "__main__":
from pylab import *
from scipy.spatial import Delaunay
print ' '
print 'Test case -- square domain with Del^2 Psi=0 and various BC'
#set up figure
figure(1)
clf()
#make points in square domain of size 1
points=[]
xvec=linspace(0.0,10.0,50)
yvec=linspace(0.0,10.0,50)
xmat,ymat=meshgrid(xvec,yvec)
#make Mask
mask=xmat<Inf
#if true, make mask irregular
makeBight=False
if makeBight:
if False: #notch
bight=logical_and(xmat<2.0,ymat>5)
mask[bight]=False
else:
#diagnonal cut,slope<1 to make jaggies
bight=(ymat>(1.5*xmat+5))
#bight=(ymat>(0.75*xmat+5))
mask[bight]=False
#now make triangulation, and convert to mesh
points,simplices,matPnts=Matrix2Delaunay(xmat,ymat,mask,addOrthog=True)
mesh=ScipyDelaunay2mesh(points,simplices)
#setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
#first, make grid variable
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
#apply boundary of Del\Psi \dot \nhat=1 on x=0 boundary
# Del\Psi \dot \nhat=0 on y=0,1 boundary
#and Psi=1 on x=1
xFace,yFace=mesh.faceCenters
faceNorm=mesh.faceNormals
figure(1)
clf()
plot(xFace[array(mesh.exteriorFaces)],yFace[array(mesh.exteriorFaces)],'co',markersize=10,mec='none')
#keep track of which boundaries have been set
boundariesNotSet=array(mesh.exteriorFaces)
#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
plot(xFace[yeq0bound],yFace[yeq0bound],'rs')
boundariesNotSet[yeq0bound]=False
#y=10.0 boundary
yeq10bound=yFace==10.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq10bound)
plot(xFace[yeq10bound],yFace[yeq10bound],'rs')
boundariesNotSet[yeq10bound]=False
#x=0 boundary
xeq0bound=xFace==0.0
#psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
psi.faceGrad.constrain(0.0*faceNorm,where=xeq0bound)
#psi.constrain(1.0,where=xeq0bound)
plot(xFace[xeq0bound],yFace[xeq0bound],'rs')
boundariesNotSet[xeq0bound]=False
#x=10 boundary
xeq10bound=xFace==10.0
psi.constrain(1.0,where=xeq10bound)
plot(xFace[xeq10bound],yFace[xeq10bound],'rs')
boundariesNotSet[xeq10bound]=False
#set remaining boundaries to zero normal gradient; this should just be bight
psi.faceGrad.constrain(0.0*faceNorm,where=boundariesNotSet)
plot(xFace[boundariesNotSet],yFace[boundariesNotSet],'ys')
#make equatiion and solve
DiffCoeff=1.0e+0
ConvecCoeff=(1.0,0.0) #flow towards x=0
#ConvecCoeff=(-1.0,0.0) #flow towards x=1
#first solve steady, then use that as intial to transient
solver=fipy.DefaultAsymmetricSolver
#start with solution to steady solution
eq=(fipy.DiffusionTerm(var=psi,coeff=DiffCoeff)+fipy.PowerLawConvectionTerm(var=psi,coeff=ConvecCoeff))
eq.solve(var=psi,solver=solver())
#now continue with transient solution. If the solution to the
#steady solution were perfect, nothing would change...
eq=(fipy.TransientTerm(var=psi)==
fipy.DiffusionTerm(var=psi,coeff=DiffCoeff)+fipy.PowerLawConvectionTerm(var=psi,coeff=ConvecCoeff))
#now loop in time and solve
dx=xvec[2]-xvec[1]
print 'for u=1 and diff=1 and dx=%f, u cfl=%f and diff cfl=%f',dx/1.0,dx**2/1.0
dt=min(dx/1.0,dx**2/1.0)*25.5
print ' and dt = ',dt
figure(2)
show()
nt=0
# while nt*dt<200.0:
while nt<10:
#and plot
if remainder(nt,1)==0:
clf()
subplot(2,1,1)
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
simplices,psi.numericValue)
colorbar()
title('solution, min=%f, max=%f at t=%3.2f'%(amin(psi.numericValue),amax(psi.numericValue),nt*dt))
axis('square')
subplot(2,1,2)
plot(array(psi.mesh.x),psi.numericValue,'r.')
draw()
pause(0.001)
print ' doing time',dt*nt,'step nt=',nt,'STD=',std(psi.numericValue)
nt+=1
eq.solve(var=psi,solver=solver(),dt=dt)
#now plot mesh
figure(1)
triplot(mesh.vertexCoords[0,:],mesh.vertexCoords[1,:],triangles=simplices)
title('simple square mesh')
axis('square')
if True:
#plot gradients at cell location
dPdx=array(psi.grad[0,:])
dPdy=array(psi.grad[1,:])
x=array(psi.mesh.x)
y=array(psi.mesh.y)
quiver(x,y,dPdx,dPdy,scale=20.0)
axis('square')
draw()
show()
#now plot solution
figure(2)
clf()
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
simplices,psi.numericValue)
colorbar()
title('solution, min=%f, max=%f'%(amin(psi.numericValue),amax(psi.numericValue)))
if False:
#plot gradients at cell location
dPdx=array(psi.grad[0,:])
dPdy=array(psi.grad[1,:])
x=array(psi.mesh.x)
y=array(psi.mesh.y)
quiver(x,y,dPdx,dPdy,scale=20.0)
axis('square')
figure(3)
clf()
plot(array(psi.mesh.x),psi.numericValue,'r.')
x=array(psi.mesh.x)
plot(x,exp(-x)+1.0,'b.',alpha=0.1)
xlabel('x')
ylabel('psi')
draw()
show()
# These modules takes the output of scipy.spatial.Delaunay's
# triangulation of data points, and returns both the FiPy mesh for the
# points, and a list of the boundary edges so that boundary conditions
# can be set for the points.
#
# This code handles the case where triangles/simplices have been
# removed from output of the Delaunay triangulation after the
# triangulation. This triangle removal can be useful to create
# interior holes in the domain, or to handle concave boundaries of the
# domain.
#
# The documentation for the routines can be found in there docstrings,
# and examples of there use can be found in __main__
from numpy import *
import fipy
from collections import defaultdict
import copy as Copy
from JMP_Make_Grids import *
if __name__ == "__main__":
from pylab import *
from scipy.spatial import Delaunay
print ' '
print 'Test case -- square domain with Del^2 Psi=0 and various BC'
#set up figure
figure(1)
clf()
#make points in square domain of size 1
points=[]
xvec=linspace(0.0,10.0,50)
yvec=linspace(0.0,10.0,50)
xmat,ymat=meshgrid(xvec,yvec)
#make Mask
mask=xmat<Inf
#if true, make mask irregular
makeBight=False
if makeBight:
if True: #notch
bight=logical_and(xmat<2.0,ymat>5)
mask[bight]=False
else:
#diagnonal cut,slope<1 to make jaggies
bight=(ymat>(1.5*xmat+5))
#bight=(ymat>(0.75*xmat+5))
mask[bight]=False
#now make triangulation, and convert to mesh
points,simplices,matPnts=Matrix2Delaunay(xmat,ymat,mask,addOrthog=True)
mesh=ScipyDelaunay2mesh(points,simplices)
#setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
#first, make grid variable
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
#apply boundary of Del\Psi \dot \nhat=1 on x=0 boundary
# Del\Psi \dot \nhat=0 on y=0,1 boundary
#and Psi=1 on x=1
xFace,yFace=mesh.faceCenters
faceNorm=mesh.faceNormals
figure(1)
clf()
plot(xFace[array(mesh.exteriorFaces)],yFace[array(mesh.exteriorFaces)],'co',markersize=10,mec='none')
#keep track of which boundaries have been set
boundariesNotSet=array(mesh.exteriorFaces)
#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
plot(xFace[yeq0bound],yFace[yeq0bound],'rs')
boundariesNotSet[yeq0bound]=False
#y=10.0 boundary
yeq10bound=yFace==10.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq10bound)
plot(xFace[yeq10bound],yFace[yeq10bound],'rs')
boundariesNotSet[yeq10bound]=False
#x=0 boundary
xeq0bound=xFace==0.0
#psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
psi.faceGrad.constrain(0.0*faceNorm,where=xeq0bound)
#psi.constrain(1.0,where=xeq0bound)
plot(xFace[xeq0bound],yFace[xeq0bound],'rs')
boundariesNotSet[xeq0bound]=False
#x=10 boundary
xeq10bound=xFace==10.0
psi.constrain(1.0,where=xeq10bound)
plot(xFace[xeq10bound],yFace[xeq10bound],'rs')
boundariesNotSet[xeq10bound]=False
#set remaining boundaries to zero normal gradient; this should just be bight
psi.faceGrad.constrain(0.0*faceNorm,where=boundariesNotSet)
plot(xFace[boundariesNotSet],yFace[boundariesNotSet],'ys')
#make equatiion and solve
DiffCoeff=1.0e+0
ConvecCoeff=(1.0,0.0) #flow towards x=0
#ConvecCoeff=(-1.0,0.0) #flow towards x=1
eq=(fipy.DiffusionTerm(var=psi,coeff=DiffCoeff)+fipy.ExponentialConvectionTerm(var=psi,coeff=ConvecCoeff))
#eq=(fipy.DiffusionTerm(var=psi,coeff=DiffCoeff)+fipy.PowerLawConvectionTerm(var=psi,coeff=ConvecCoeff))
#eq=(fipy.DiffusionTerm(var=psi,coeff=DiffCoeff)+fipy.CentralDifferenceConvectionTerm(var=psi,coeff=ConvecCoeff))
#eq=(fipy.DiffusionTerm(var=psi,coeff=DiffCoeff))
eq.solve(var=psi,solver=fipy.LinearLUSolver())
#now plot mesh
triplot(mesh.vertexCoords[0,:],mesh.vertexCoords[1,:],triangles=simplices)
title('simple square mesh')
axis('square')
if True:
#plot gradients at cell location
dPdx=array(psi.grad[0,:])
dPdy=array(psi.grad[1,:])
x=array(psi.mesh.x)
y=array(psi.mesh.y)
quiver(x,y,dPdx,dPdy,scale=20.0)
axis('square')
draw()
show()
#now plot solution
figure(2)
clf()
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
simplices,psi.numericValue)
colorbar()
title('solution, min=%f, max=%f'%(amin(psi.numericValue),amax(psi.numericValue)))
if False:
#plot gradients at cell location
dPdx=array(psi.grad[0,:])
dPdy=array(psi.grad[1,:])
x=array(psi.mesh.x)
y=array(psi.mesh.y)
quiver(x,y,dPdx,dPdy,scale=20.0)
axis('square')
figure(3)
clf()
plot(array(psi.mesh.x),psi.numericValue,'r.',label='Numerical Solution')
x=array(psi.mesh.x)
plot(x,exp(-x)+1.0,'b.',alpha=1.,label='1+exp(-x)')
legend()
xlabel('x')
ylabel('psi')
draw()
show()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
