Hi Raymond, I was just corresponding with James Pringle off the list about the same issue for a Delaunay triangulated mesh. Here is my explanation to him.
I think this is the non-orthogonality issue. I should have realized this much sooner. The only place we seem to reference this issue is in www.ctcms.nist.gov/fipy/documentation/numerical/discret.html where it is stated "this estimate relies on the orthogonality of the mesh, and becomes increasingly inaccurate as the non-orthogonality increases. Correction terms have been derived to improve this error but are not currently included in FiPy [13]." See this for example, https://www.researchgate.net/publication/242349760_Finite_volume_method_for_the_solution_of_flow_on_distorted_meshes In the test mesh that you show, some of the triangles are highly non-orthogonal to each other. This results in errors. You can view the quantity of non-orthogonality in the mesh using nonorth_var = fipy.CellVariable(mesh=mesh, value=mesh._nonOrthogonality) fipy.Viewer(nonorth_var).plot() The following code uses triangles, but with a perfectly orthogonality and, thus, gives perfect results. ~~~~ import fipy as fp nx = 50 dx = 1. / nx mesh = fp.Tri2D(nx=nx, ny=nx, dx=dx, dy=dx) psi = fp.CellVariable(mesh=mesh) psi.constrain(1.0, where=mesh.facesLeft) psi.faceGrad.constrain(1.0, where=mesh.facesRight) eq = fp.DiffusionTerm().solve(psi) fp.Viewer(psi - mesh.cellCenters[0]).plot() nonorth_var = fp.CellVariable(mesh=mesh, value=mesh._nonOrthogonality) fp.Viewer(nonorth_var).plot() raw_input('stopped') ~~~~ As far as your Delaunay triangulation is concerned, finding a way to smooth out / reduce the non-orthogonality would be one way to improve the accuracy. The FiPy docs definitely need to be more explicit about this issue. Also, actually addressing the orthogonality issue for diffusion terms at least would be good. Cheers, Daniel On Wed, Jul 20, 2016 at 1:30 PM, Raymond Smith <smit...@mit.edu> wrote: > Hi, FiPy. > > I was looking over the diffusion term documentation, > http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#diffusion-term > and I was wondering, do we lose second order spatial accuracy as soon as we > introduce any non-uniform spacing (anywhere) into our mesh? I think the > equation right after (3) for the normal component of the flux is only second > order if the face is half-way between cell centers. If this does lead to > loss of second order accuracy, is there a standard way to retain 2nd order > accuracy for non-uniform meshes? > > I was playing around with this question here: > https://gist.github.com/raybsmith/e57f6f4739e24ff9c97039ad573a3621 > with output attached, and I couldn't explain why I got the trends I saw. > The goal was to look at convergence -- using various meshes -- of a simple > diffusion equation with a solution both analytical and non-trivial, so I > picked a case in which the transport coefficient varies with position such > that the solution variable is an arcsinh(x). I used three different styles > of mesh spacing: > * When I use a uniform mesh, I see second order convergence, as I'd expect. > * When I use a non-uniform mesh with three segments and different dx in each > segment, I still see 2nd order convergence. In my experience, even having a > single mesh point with 1st order accuracy can drop the overall accuracy of > the solution, but I'm not seeing that here. > * When I use a mesh with exponentially decreasing dx (dx_i = 0.96^i * dx0), > I see 0.5-order convergence. > > I can't really explain either of the non-uniform mesh cases, and was curious > if anyone here had some insight. > > Thanks, > Ray > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Daniel Wheeler _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
