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 <[email protected]> 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 > (http://www.ctcms.nist.gov/fipy/fipy/generated/fipy.variables.html#fipy.variables.cellVariable.CellVariable.faceGrad) > 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) > <[email protected]> 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 <[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] > > 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 > [email protected] > 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 > [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 ]
