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