Re: scipy's Delaunay output to fipy mesh2D object

[James Pringle]

Dear all --

A warning about using the above code. Since the choice of triangles
formed by a Delaunay triangulation can be somewhat arbitrary in the case of
a somewhat regular distribution of points, it will choose some
triangulations that lead to severe orthogonality errors in fipy, even if
that was not necessary for the points it was given.

   So use with caution and monitor non-orthogonality with the
mesh._nonOrthogonal property. For my mostly regular grid I am resorting to
triangulation-by-hand.

Cheers,
Jamie Pringle
University of New Hampshire

On Fri, Jun 17, 2016 at 10:14 AM, James Pringle  wrote:

> Thanks --
>
>Once I have pounded on this a bit more, I will write the function. It
> should be straightforward. I concentrated on making it simple and easy to
> understand; since it is run only to make the grid, efficiency was not a
> great issue for me.
>
> Jamie
>
> On Fri, Jun 17, 2016 at 10:03 AM, Daniel Wheeler <
> daniel.wheel...@gmail.com> wrote:
>
>> James, this is awesome, thanks for posting. It would be a very good
>> idea to have a DelaunayMesh in FiPy that used Scipy for the
>> triangulation as it would give another way to make and test
>> unstructured meshes without the Gmsh dependency. Something like,
>>
>> mesh = DelaunayMesh(points_in_2D)
>>
>> In the long run, it would be good to change your code so that it
>> doesn'y loop in Python. I suspect that your code is quite inefficient
>> due to these loops. However, you have something that works, which is
>> great and it is only a one time cost at the beginning of the
>> simulation. Is efficiency an issue for you with this?
>>
>> On Thu, Jun 16, 2016 at 4:29 PM, James Pringle  wrote:
>> > For the use of others --
>> >
>> >  Sometimes when dealing with arbitrary grids, it is useful to use
>> the
>> > output of routines other than Gmsh to create the fipy mesh. In my case,
>> I
>> > need to create a complex grid from ocean bathymetry, and I ended up
>> using
>> > scipy.spatial.Delaunay to do so. With help from Daniel and Jonathan, I
>> > created a code to create a mesh2D object from the output of Delaunay.
>> It is
>> > written to emphasize clarity over speed. I hope others find it useful.
>> >
>> > James Pringle
>> >
>> > #==
>> > #This code creates a fipy grid from the output of a Delaunay
>> > #triangulation. The code is written to prioratize clarity over speed.
>> >
>> > from pylab import *
>> > from numpy import *
>> > import fipy
>> > from scipy.spatial import Delaunay
>> >
>> >
>> > #make a simple example set of points to triangulate with delaunay
>> > points=array([[ 0.,  0.],
>> >[ 0.,  1.],
>> >[ 1.,  0.],
>> >[ 1.,  1.]])
>> > tri=Delaunay(points)
>> >
>> > #when this code is run as is, it should produce
>> > #  faceVertexIDs=
>> > #   [[3, 1, 0, 2, 0
>> > #[1, 0, 3, 3, 2]]
>> > #and
>> > #  cellFaceIds=
>> > #[[0, 3],
>> > # [1, 2],
>> > # [2, 4]]
>> >
>> >
>> > #now 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 tri.points in
>> > #content, but reshaped.
>> > print 'Making vertexCoords'
>> > vertexCoords=tri.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(tri.simplices.shape[0]):
>> > for nface in xrange(3):
>> > if nface==0:
>> > face=tri.simplices[ntri,[0,1]]
>> > elif nface==1:
>> > face=tri.simplices[ntri,[1,2]]
>> > else:
>> > face=tri.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 

Re: Accuracy of DiffusionTerm for non-uniform mesh

[Raymond Smith]

Thanks.

And no, I'm not sure about the normalization for grid spacing. I very well
could have calculated the error incorrectly. I just reported the root mean
square error of the points and didn't weight by volume or anything like
that. I can change that, but I'm not sure of the correct approach.

On Wed, Jul 20, 2016 at 4:25 PM, Daniel Wheeler 
wrote:

> On Wed, Jul 20, 2016 at 1:30 PM, Raymond Smith  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?
>
> This is a different issue than the non-orthogonality issue, my mistake
> in the previous reply.
>
> > 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.
>
> That's strange. Are you sure that all the normalization for grid
> spacing is correct when calculation the norms in that last case?
>
> > I can't really explain either of the non-uniform mesh cases, and was
> curious
> > if anyone here had some insight.
>
> I don't have any immediate insight, but certainly needs to addressed.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Accuracy of DiffusionTerm for non-uniform mesh

[Daniel Wheeler]

On Wed, Jul 20, 2016 at 1:30 PM, Raymond Smith  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?

This is a different issue than the non-orthogonality issue, my mistake
in the previous reply.

> 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.

That's strange. Are you sure that all the normalization for grid
spacing is correct when calculation the norms in that last case?

> I can't really explain either of the non-uniform mesh cases, and was curious
> if anyone here had some insight.

I don't have any immediate insight, but certainly needs to addressed.

-- 
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Improving non-orthogonality issue in FiPy

[Daniel Wheeler]

FiPy has issues with non-orthogonal meshes (i.e. when the vector
between the cell center isn't parallel to the face normal). We did
make an attempt at fixing this, which resulted in two diffusion terms,


https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermCorrection.py

and


https://github.com/usnistgov/fipy/blob/develop/fipy/terms/diffusionTermNoCorrection.py

As I recall, the "corrected" diffusion term was quite unstable for
moderate to high non-orthogonality.

I'd like to address this issue again and I was wondering if anyone had
opinions regarding the latest / best / most tractable approach for
dealing with non-orthogonality in CC-FV.  I haven't followed the
literature on this for quite some time.

-- 
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Accuracy of DiffusionTerm for non-uniform mesh

[Daniel Wheeler]

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  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 ]

Accuracy of DiffusionTerm for non-uniform mesh

[Raymond Smith]

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_convergence.pdf
Description: Adobe PDF document
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]