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