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 fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]