Ok, Then there is also a subtle issue with psi.faceGrad or something else.
As above, I am solving the problem
eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
With the BC's
\Del\Psi \dot \nhat = -1 at x=0
\Del\Psi \dot \nhat =0 at y=0 and y=1
Psi =1 at x=1
With the code
#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
#y=1.0 boundary
yeq1bound=yFace==1.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
#x=0 boundary
xeq0bound=xFace==0.0
psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
#x=1 boundary
xeq1bound=xFace==1.0
psi.constrain(0.0,where=xeq1bound)
The solution on the cell variables matches the analytical solution of
psi=x, but when I get the gradient with psi.faceGrad it is incorrect on the
y=0 and y=1 boundary (it is correct everywhere else, including the x
boundaries). On the y boundaries, psi.faceGrad returns (0,0) at each point,
instead of the correct (1,0). This is not such a big deal, but should be
documented.
Code illustrating the error is attached.
Thanks,
Jamie
On Mon, Jul 18, 2016 at 3:17 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:
> They return \nabla\psi, the gradient, not the projection onto the normal
> (or anything else).
>
> > On Jul 18, 2016, at 1:47 PM, James Pringle <jprin...@unh.edu> wrote:
> >
> > Dear Jonathan --
> >
> > I have to admit I picked faceGradAverage over faceGrad randomly. One
> comment on my test code -- it shuffles the order of the points that
> Delaunay uses to make the mesh each run. I did this as a debugging step of
> my Delaunay code. If the issue is the order of the triangle points, you
> should get different answers each time, and indeed you do. Just search
> for "shuffle" in my test code.
> >
> > Also, the documentation for faceGrad and faceGradAverage (
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy_fipy_generated_fipy.variables.html-23fipy.variables.cellVariable.CellVariable.faceGrad&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=vHBVDYOH9hZ7yTR_To6QH5i-5bO4QbHKmXpEShgMJ54&e=
> ) is unclear. Do these routines return \Del\psi, or \Del\psi\dot\nhat?
> I.e. do they return the gradient or the gradient projected onto the unit
> normal to the face?
> >
> > If it is the former, there is also a problem in faceGrad.
> >
> > Thanks!
> > Jamie
> >
> > On Mon, Jul 18, 2016 at 1:23 PM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
> > James -
> >
> > I didn't even know we had a .faceGradAverage!
> >
> > The issue turns out to be with .grad (.faceGradAverage is calculated
> from .grad and .grad is calculated from .faceValue, whereas .faceGrad is
> calculated from .value (at cell centers)). Neither .grad nor
> .faceGradAverage is used in calculating the solution to the Laplace
> problem, which is why the solution looks OK.
> >
> > My guess is that some of the triangles you send to Mesh2D are
> "upside-down". We don't stipulate an ordering and that upside-downess
> should impact the calculation of the diffusion stencil, so our best guess
> is that we correct for the vertex-face-cell ordering elsewhere and just
> need to apply the same correction to .grad. We'll try to sort this out in
> the next few days.
> >
> > - Jon
> >
> > > On Jul 15, 2016, at 3:20 PM, James Pringle <jprin...@unh.edu> wrote:
> > >
> > > I am trying to debug some boundary conditions in a complex problem. I
> have run into a problem with faceGradAverage or my understanding of it. It
> is illustrated in a simple problem
> > >
> > > \Del^2 Psi =0 on a unit square
> > > normal derivative = 1 on x=0
> > > normal derivative = 0 on y=0
> > > normal derivative = 0 on y=1
> > > Psi =1 on x=1
> > >
> > > The analytical solution to this problem is
> > >
> > > Psi = x
> > >
> > > The code is in the attached file; skip until "if __name__ ==
> "__main__":" , the earlier code converts a Delaunay triangulation to a fipy
> mesh. I am doing this to match my more complex case.
> > >
> > > Anyway, when the code is run, I get the analytical solution, and so
> gradient \Del\Psi should be (1,0) everywhere. It is not. Note also that
> since \Del\Psi is spatially uniform, how the gradient is averaged in space
> should not make a difference.
> > >
> > > How am I confused?
> > >
> > > So that you don't have to open the attachment, the code that solves
> the problem and prints the face gradients is included as text here:
> > >
> > > #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
> (equivalent to d Psi/dx=1)
> > > # Del\Psi \dot \nhat=0 on y=0,1 boundary
> > > #and Psi=1 on x=1
> > >
> > > xFace,yFace=mesh.faceCenters
> > > faceNorm=mesh.faceNormals
> > >
> > > #y=0 boundary
> > > yeq0bound=yFace==0.0
> > > psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
> > >
> > > #y=1.0 boundary
> > > yeq1bound=yFace==1.0
> > > psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
> > >
> > > #x=0 boundary
> > > xeq0bound=xFace==0.0
> > > psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
> > >
> > > #x=1 boundary
> > > xeq1bound=xFace==1.0
> > > psi.constrain(0.0,where=xeq1bound)
> > >
> > >
> > > #make equatiion and solve
> > > eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
> > > eq.solve(var=psi)
> > >
> > > #now plot solution
> > > subplot(2,1,2)
> > > tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
> > > tri.simplices,psi.numericValue)
> > > colorbar()
> > > title('solution')
> > >
> > > #now print faceGradients
> > > print 'Face Gradients on x=0, should be
> 1,0',psi.faceGradAverage[:,xeq0bound]
> > > print ' '
> > > print 'Face Gradients on y=0, should be
> 1,0',psi.faceGradAverage[:,yeq0bound]
> > >
> > > Thanks,
> > > Jamie
> > > <SquareGrid_Clean.py>_______________________________________________
> > > fipy mailing list
> > > fipy@nist.gov
> > >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=10-koYdkMl_KoP5wN4ux7A9dDrjpgu_k8oMKcMpXxAk&s=msbLTmmKWXFW9nc95yGwb31rXfgwljQh5qj51_e-6Eg&e=
> > > [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=10-koYdkMl_KoP5wN4ux7A9dDrjpgu_k8oMKcMpXxAk&s=VR9IvO4aq-1MzIiDS7oISwOtQO-TgcN8rowVOl5J_wE&e=
> ]
> >
> >
> > _______________________________________________
> > fipy mailing list
> > fipy@nist.gov
> >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=10-koYdkMl_KoP5wN4ux7A9dDrjpgu_k8oMKcMpXxAk&s=msbLTmmKWXFW9nc95yGwb31rXfgwljQh5qj51_e-6Eg&e=
> > [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=10-koYdkMl_KoP5wN4ux7A9dDrjpgu_k8oMKcMpXxAk&s=VR9IvO4aq-1MzIiDS7oISwOtQO-TgcN8rowVOl5J_wE&e=
> ]
> >
> > _______________________________________________
> > fipy mailing list
> > fipy@nist.gov
> >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=m5I6LErPy5RM6awOYyB3o_2cvOQ7P8bR5_7OyR2-6h8&e=
> > [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=sXQupvT3I1wAAqKu1Jm3gkRW-puZq01lo-1q61L3vd4&e=
> ]
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=m5I6LErPy5RM6awOYyB3o_2cvOQ7P8bR5_7OyR2-6h8&e=
> [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=sXQupvT3I1wAAqKu1Jm3gkRW-puZq01lo-1q61L3vd4&e=
> ]
>
# 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
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
borderList: a list of lists
Each list in borderList is a list of the indices of the faces
as indexed in mesh.cellFaceIDs. Each list is for a single
contiguous border of the domain. These will be used to set the
boundary conditions on the mesh.
'''
#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)
#now make borderList. We cannot assume the neighbors property of
#the Delaunay triangulation is valid because some triangles may
#have been removed.
#first make vertexPathList, a list of all points that are on
#borders
#first, redo neighbors. Create a defaultdict that when called
#with a new key returns an empty list
whichTriangles=defaultdict(list);
#the entry into the dictionary are the edges of the triangle, sorted
#so smaller vertex is first
print ' Making new neighbor data'
for n in range(simplices.shape[0]):
edge1=tuple(sort(simplices[n,[0,1]]))
edge2=tuple(sort(simplices[n,[1,2]]))
edge3=tuple(sort(simplices[n,[2,0]]))
whichTriangles[edge1].append(n)
whichTriangles[edge2].append(n)
whichTriangles[edge3].append(n)
if debug:
triplot(points[:,0],points[:,1],triangles=simplices)
#to debug, plot over all triangles and plot there number in the centers
for n in range(simplices.shape[0]):
xave=0.0
yave=0.0
for ne in range(3):
xave=xave+points[simplices[n,ne],0]/3.0
yave=yave+points[simplices[n,ne],1]/3.0
text(xave,yave,str(n),fontsize=6)
#to debug, iterate over all edges in whichTriangle, and plot number of neighbors
for edge in whichTriangles.keys():
xedge=0.5*(points[edge[0],0]+points[edge[1],0])
yedge=0.5*(points[edge[0],1]+points[edge[1],1])
if len(whichTriangles[edge])==1:
plot(xedge,yedge,'r*',mec='none')
text(points[edge[0],0],points[edge[0],1],str(edge[0]),color='r',fontsize=9)
text(points[edge[1],0],points[edge[1],1],str(edge[1]),color='r',fontsize=9)
elif len(whichTriangles[edge])==2:
plot(xedge,yedge,'g*',mec='none')
else:
assert False, 'wrong number of neighbors'
#now, for each edge/face, see if it has any neighbors
neighbors=simplices*0
for n in range(simplices.shape[0]):
edge1=tuple(sort(simplices[n,[0,1]]))
edge2=tuple(sort(simplices[n,[1,2]]))
edge3=tuple(sort(simplices[n,[2,0]]))
edges=[edge1,edge2,edge3]
for ne in range(3):
edge=edges[ne]
if len(whichTriangles[edge])==1:
neighbors[n,ne]=-1
elif len(whichTriangles[edge])==2:
if whichTriangles[edge][0]==n:
neighbors[n,ne]=whichTriangles[edge][1]
else:
neighbors[n,ne]=whichTriangles[edge][0]
assert whichTriangles[edge][1]==n,'Umm, how can this be true?!?!'
else:
assert False, 'I am confused, simplex has other than 0 or 1 neighbor...'
#identify border edges/face; they are the ones where
#neighbor=-1. Collect edgePoints, a list of tuples that define
#points that connect along the edges of the domain
edgePoints=[]
for n in range(simplices.shape[0]):
edge1=sort(simplices[n,[0,1]])
edge2=sort(simplices[n,[1,2]])
edge3=sort(simplices[n,[2,0]])
if neighbors[n,0]==-1:
edgePoints.append((edge1[0],edge1[1]))
if debug:
plot(points[[edge1[0],edge1[1]],0],
points[[edge1[0],edge1[1]],1],'r-',linewidth=4,alpha=0.2)
if neighbors[n,1]==-1:
edgePoints.append((edge2[0],edge2[1]))
if debug:
plot(points[[edge2[0],edge2[1]],0],
points[[edge2[0],edge2[1]],1],'r-',linewidth=4,alpha=0.2)
if neighbors[n,2]==-1:
edgePoints.append((edge3[0],edge3[1]))
if debug:
plot(points[[edge3[0],edge3[1]],0],
points[[edge3[0],edge3[1]],1],'r-',linewidth=4,alpha=0.2)
#ok, now make vertexPathList. Make a dictionary of paths so
#paths(npoint) returns list of points the path is connected to.
print ' Now make vertexPathList'
paths=defaultdict(list)
for np in xrange(len(edgePoints)):
edge=edgePoints[np]
paths[edge[0]].append(edge[1])
paths[edge[1]].append(edge[0])
#now follow each path to find a cycle; store each cycle in vertexPathList,
#which is a list of lists. Be sure to remove each path after taking it.
vertexPathList=[]
keys=paths.keys()
for key in keys:
edges=paths[key]
if len(edges)==0:
assert True,'nothing to do here'
else:
thispath=[key] #start of path
nhere=paths[key].pop()
nfrom=key+0
cols='rgbcmyk'
while (nhere != key):
if False:
plot([points[nhere,0],points[nfrom,0]]
,[points[nhere,1],points[nfrom,1]],cols[remainder(key,7)]+'*-',
mec='none',markersize=12,linewidth=2)
#remove path whence we came
paths[nhere].remove(nfrom)
#add this point to path
thispath.append(nhere)
#now move to next point, also remove that path:
nfrom=nhere+0
nhere=paths[nhere].pop()
vertexPathList.append(thispath)
#vertexPathList contains a list of lists of vertex points that
#surround all islands and (at item nOcean) the model domain
#boundary. All of paths are closed. I need to convert this into a
#list of faces, as defined in faceMap and faceVertexIDs. I also
#need to make sure to include the face that connects the first and
#last element of the list.
#remember faceMap is a dictionary with a key made up of a sorted tuple of vertex
#indices which maps from a face to its location in faceVertexIDs
#first, make a copy of vertexPathList, since all lists will be the same size
borderList=Copy.deepcopy(vertexPathList)
#now fill fipy BoundaryList with good data
for ncycle in xrange(len(vertexPathList)):
for npoint in xrange(len(vertexPathList[ncycle])):
nOtherPoint=remainder(npoint+1,len(vertexPathList[ncycle]))
faceKey=tuple(sort((vertexPathList[ncycle][npoint],
vertexPathList[ncycle][nOtherPoint])))
borderList[ncycle][npoint]=faceMap[faceKey]
return (mesh,borderList)
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,1.0,30)
xmat,ymat=meshgrid(xvec,xvec)
points=zip(xmat.flatten(),ymat.flatten())
#to debug shuffle points
shuffle(points)
#now make triangulation, and convert to mesh
tri=Delaunay(points)
mesh,borderList=ScipyDelaunay2mesh(tri.points,tri.simplices)
#now plot mesh
figure(1)
clf()
subplot(2,1,1)
triplot(mesh.vertexCoords[0,:],mesh.vertexCoords[1,:],triangles=tri.simplices)
title('simple square mesh')
draw()
show()
#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
#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
#y=1.0 boundary
yeq1bound=yFace==1.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
#x=0 boundary
xeq0bound=xFace==0.0
psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
#x=1 boundary
xeq1bound=xFace==1.0
psi.constrain(0.0,where=xeq1bound)
#make equatiion and solve
eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
eq.solve(var=psi)
#now plot solution
subplot(2,1,2)
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
tri.simplices,psi.numericValue)
colorbar()
title('solution')
#now print faceGradients
print 'Face Gradients on x=0, should be 1,0',psi.faceGrad[:,xeq0bound]
print ' '
print 'Face Gradients on y=0, should be 1,0',psi.faceGrad[:,yeq0bound]
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
