Buggy behaviour on macosx using enthought?

[Jon Alm Eriksen]

Hi there,

Thanks for providing this code, this is exactly what I need for 2-phase flow simulation. I'm not an expert on front tracking, so I might have this completely wrong. In any case, I need to approximate the curvature of a zero-contour. I calculate the distance function (order=2) to solve the eikonal equation given the zero-contour. Then I use a 9pt laplacian stencil to approximate the curvature. Although the distance function seems to give reasonable results, I get that |grad distance| evaluated by central differences at certain points deviate significantly from 1. I believe that this error also makes the curvature approximation completely off (see contour plot below).

Here is a minimal script which illustrates my problem:

http://folk.uio.no/jonaerik/files/minimal.py

the output:

evaluating Q = ||grad phi| - 1| near zero-contour:

max(Q) = 0.145076690953

mean(Q) = 0.0127694905239

grid spaceing = 0.00500500500501

grid shape = (1000, 1000)

the countor plot:

http://folk.uio.no/jonaerik/files/example.png

This might be related to the following test error, given by 

python -c "import skfmm; skfmm.test()" > testerror.txt

http://folk.uio.no/jonaerik/files/testerror.txt

I use the latest version of Enthought Python Distribution (EPD) on a Macbook Pro using OSX 10.8.4

Is the output of my script related to the test error, and if so, is there any way of fixing the behaviour? Any help would be highly appreciated.

best

Jon Alm Eriksen

[Jason Furtney]

Jon,

Thanks for the email with this information. I think the magnitude of the gradient is noisy in this example because the initial zero level set of phi has sharp angles. If you define the initial contour as 

phi = 0.5*X**2 + Y**2 -1

the magnitude of the gradient is more uniform.  This makes your curvature calculation better but there are still some problems. The curvature kappa is usually written as div (grad phi/|grad phi|) so maybe normalizing the gradient will give better results? I have never worked with front curvature so I am not 100% sure. I know the Sethian book has a chapter about front curvature. 

Another issue is that the starting narrow band initialization is first order, this could be improved. This would probably help in this case because all the action is near the initial front.

The error in the fast marching method seems to propagate along the grid diagonals. There is a multi-stencil method which improves this. http://www.cvip.uofl.edu/wwwcvip/research/publications/Pub_Pdf/2007_2/Sabry_Farag_TPAMI_2007.pdf but it is not implemented in this module. 

About the test failures, these are some corner cases we are working on improving, as far as I know the module is working as intended. 

Have a look at lsmlib and the (pylsmlib wrappers) as it may have better front curvature functionality. 

I hope that helps, I will give this some more thought and get back to you if I come up with anything else. If you get some nice front curvature plots it would be great to include them as an example in the scikit-fmm documentation. 

Thanks and let me know how you get on,

Jason

Jon Alm Eriksen

--
You received this message because you are subscribed to the Google Groups "scikit-fmm" group.
To unsubscribe from this group and stop receiving emails from it, send an email to scikit-fmm+...@googlegroups.com.
To post to this group, send email to sciki...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/816d3c83-61cd-4094-b22d-1aa974cd60e3%40googlegroups.com?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

--
--
  Jason K. Furtney
  Itasca Consulting Group
  111 3rd Ave. South, Suite 450
  Minneapolis, MN 55401 USA
    (612) 371-4711

[Jon Alm Eriksen]

Thank you so much for the answer Jason, here comes a couple of
followup questions/comments.

> ... I think the magnitude of the

> gradient is noisy in this example because the initial zero level set of phi
> has sharp angles. If you define the initial contour as
>
> phi = 0.5*X**2 + Y**2 -1

My understanding of the FMM was that the distance function was
uniquely determined by the algorithm, given the identification of
nodes inside and outside a boundary. At least, that is how I
understood the chapter on FMM in Sethian's book. If this is not the
case, how do you use the initial values to determine phi, are you
using some kind of front reconstruction scheme to find the position of
the front between inside and outside nodes?

> .... The curvature kappa is

> usually written as div (grad phi/|grad phi|) so maybe normalizing the
> gradient will give better results?

This is the general expression for the curvature; if you evaluate the
expression at a point, it will correspond to the curvature of the
contour of that point. If phi is a solution of the eikonal equation,
we have by definition that |grad phi| = 1, and the above expression
will simplify to the laplacian of phi. I don't think that using the
general expression would improve the results. The unwanted behaviour
is related to that |grad phi| isn't close enough to 1, and I see no
reason why it would reduce the error by introducing an approximation
of |grad phi| in the expression.

Another thing I don't understand is that the code doesn't seem to
converge to the analytical solution as I increase the grid
size/decrease the lattice spacing. Rather, the error takes a finite
value (e.g. max(Q) in my code example). This might be a property of
the scheme. To be able to extract the laplacian of phi, I guess the
numerical scheme has to be accurate to a higher order. This article
http://arxiv.org/abs/1111.6903 suggest an approach, but I don't know
how hard it would be to implement. (btw what exactly is the 'order'
parameter input referring to?)

Also thanks for the suggestions on alternative implementations. I had
some problems getting lsmlib to work properly. I just installed FiPy,
and will give it a try.

If I get the skfmm package to work properly, I would be more than
happy about writing a tutorial example on curvature.


best
jon

[Daniel Wheeler]

On Mon, May 20, 2013 at 7:14 PM, Jon Alm Eriksen <jon.alm...@gmail.com> wrote:

Also thanks for the suggestions on alternative implementations. I had
some problems getting lsmlib to work properly. I just installed FiPy,
and will give it a try.

PyLSMLIB should (hopefully) give exactly the same results as Skfmm as they have identical test suites now.

--
Daniel Wheeler

[Jason Furtney]

Jon and Daniel,

Thanks for the replies.

The function distanceMarcher::initalizeFrozen() in
skfmm/distance_marcher.cpp finds values for the front distance of
points adjacent to the zero contour by using bi-linear interpolation.

I think the problem is this initial part is first order, so when you
take the Laplacian of the resulting distance function the error is
amplified. I think this initializeFrozen() function can be made second
order and this should improve the curvature results. The method should
converge as the number of points goes up, if not there is a problem.

I will have a look at the paper, at first glance it looks reasonable.
They use bi-cubic interpolation in the initialization (section 2.2).

The order optional argument is the order of the spatial derivatives
used in the point updates.

Thanks,
Jason

> --
> You received this message because you are subscribed to the Google Groups "scikit-fmm" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to scikit-fmm+...@googlegroups.com.
> To post to this group, send email to sciki...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/CAGaaLYOKTZuC6e9mZex%3DjE2yL3s3AOsNsW0BnxA6NM%3D6V-iQmA%40mail.gmail.com?hl=en.

[Jon Alm Eriksen]

Jason and Daniel,

Thanks for your replies. The convergence issue seems only to affect
nodes close to the boundary (phi=0). When I choose other nodes exept
for the 'band' nodes in my code example, max(Q) goes to zero with
increasing resolution, as expected.

Let me know if you implement second order accuracy in the
initialization (bi-cubic interpolation). That should do the magic for
evaluating curvature at the front.

Also Jason, you might be right about using div(grad phi/|grad phi|)
for the curvature, even in the case of a signed distance function. It
seems like that is the preffered choice in the literature and it
might help the accuracy.

best
jon

> You received this message because you are subscribed to a topic in the Google Groups "scikit-fmm" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/scikit-fmm/ZqYjm1ii1iI/unsubscribe?hl=en.
> To unsubscribe from this group and all its topics, send an email to scikit-fmm+...@googlegroups.com.

> To post to this group, send email to sciki...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/CAOU2ZSX2uC04K%3DasN%2BQY3BAT3HqMBuHUVgghWEL6GsydSA3RRw%40mail.gmail.com?hl=en.