Seemingly erroneous 0 in travel time

[milisims]

Take the following phi and speed arrays:

import numpy as np
import skfmm
phi = np.array([[ 1,  1, -1, -1, -1, -1,  1],
                [ 3,  1, -1, -1, -1, -1,  1],
                [ 3,  3,  1,  1,  1,  1,  3],
                [ 3,  3,  3,  3,  3,  3,  3],
                [ 3,  3,  3,  3,  3,  3,  3]])
speed = np.array([[ 1.   ,  1.   ,  0.057,  0.128,  0.037,  0.039,  1.   ],
                  [ 1.   ,  1.   ,  0.199,  0.408,  0.997,  0.688,  1.   ],
                  [ 1.   ,  1.   ,  1.   ,  1.   ,  1.   ,  1.   ,  1.   ],
                  [ 1.   ,  1.   ,  1.   ,  1.   ,  1.   ,  1.   ,  1.   ],
                  [ 1.   ,  1.   ,  1.   ,  1.   ,  1.   ,  1.   ,  1.   ]])

The returned value for the fourth element in the first row when executing skfmm.travel_time(phi, speed) is given as zero, even though it is 'surrounded' by inaccessible values with speeds less than 1.

My output is:

In : skfmm.travel_time(phi, speed)
Out: 
array([[  1.5  ,   0.5  ,   8.772,   0.   ,  11.557,  12.821,   0.5  ],
       [  1.5  ,   0.5  ,   1.777,   1.225,   0.502,   0.514,   0.5  ],
       [  1.748,   0.971,   0.5  ,   0.5  ,   0.5  ,   0.5  ,   0.971],
       [  2.301,   1.748,   1.5  ,   1.5  ,   1.5  ,   1.5  ,   1.748],
       [  3.055,   2.655,   2.5  ,   2.5  ,   2.5  ,   2.5  ,   2.655]])

In : skfmm.distance(phi, 1)
Out : 
array([[ 1.499,  0.5  , -0.5  , -1.305, -1.305, -0.5  ,  0.5  ],
       [ 1.451,  0.5  , -0.354, -0.5  , -0.5  , -0.354,  0.5  ],
       [ 1.762,  0.971,  0.5  ,  0.5  ,  0.5  ,  0.5  ,  0.971],
       [ 2.337,  1.762,  1.451,  1.499,  1.499,  1.451,  1.762],
       [ 3.099,  2.673,  2.435,  2.498,  2.498,  2.435,  2.673]])

[jkfurtney]

Hi and thanks for getting in touch with this bug, I have not seen this one before.

It looks like this issue is similar to #15. In this case, the discriminant of the quadratic solved to update a point is negative and this condition is not checked for properly. I added a test for this in 5a1e1db and an exception is thrown. I will look in more detail at why this is happening this coming week.

[jkfurtney]

It looks like this is related to the second order point update. Doing the calculation in first order seems to work:

print skfmm.travel_time(phi, speed, order=1)
array([[ 1.5, 0.5, 8.77192982, 9.03360658, 23.39632103, 12.82051282, 0.5 ],
       [ 1.5, 0.5, 1.7766502, 1.2254902, 0.50150451, 0.51388574, 0.5 ],
       [ 2.04532893, 1.20710678, 0.5, 0.5, 0.5, 0.5, 1.20710678],
       [ 2.75243571, 2.04532893, 1.5, 1.5, 1.5, 1.5, 2.04532893],
       [ 3.54804305, 2.94223041, 2.5, 2.5, 2.5, 2.5, 2.94223041]])

I will keep digging.

[jkfurtney]

OK, a little more information. It looks like this bug goes back to the first release of this package. At least now we know it was not introduced recently.

[jkfurtney]

Here is a simpler version which shows the problem

import numpy as np
import skfmm
phi = np.array([[1, -1, -1],
                [0,  0, -1],
                [0,  0,  1]])

speed = np.array([[ 1, 0.01, 0.1],
                  [ 0, 0,    0.1],
                  [ 0, 0,    1. ]])
phi = np.ma.MaskedArray(phi,phi==0)
print skfmm.travel_time(phi,speed)

[jkfurtney]

Below, when value = 0.0660411041104, we get an answer. When value = 0.0660408040804, we get the negative discriminant error message. The arrival time in cell [0,2] goes asymptotically to zero between these values (which is non-physical).

phi = np.array([[1, -1, -1],
                [0,  0, -1],
                [0,  0,  1]])

speed = np.array([[ 1, value, 0.1],
                  [ 0, 0,    1. ],
                  [ 0, 0,    1. ]])

phi = np.ma.MaskedArray(phi,phi==0)
print skfmm.travel_time(phi,speed)

This only happens when two orthogonal parts of two points of the second-order stencil cross the zero level-set in the travel_time function. I am thinking it is not physically correct to construct the gradient this way. Maybe there is a different way to construct a second-order gradient or maybe we need to revert to first-order in these cases?

[jkfurtney]

The first-order method has the same problem, just at different values. Around value = 0.03405940594059409 in this case.

[jkfurtney]

The zero level set of phi going through a region with a high contrast in speed seems to cause this problem. It looks like this problem is happening when points adjacent to the initially frozen points are updated for the first time, before the main marching algorithm starts. In this simple example,

import numpy as np
import skfmm
phi = np.array([[1, -1, -1],
                [0,  0, -1],
                [0,  0,  1]])

speed = np.array([[ 1, value, 1],
                  [ 0, 0,     1],
                  [ 0, 0,     1]])

phi = np.ma.MaskedArray(phi,phi==0)
tt =  skfmm.travel_time(phi,speed,order=order)

as value decreases from 1 the arrival time at point (0,2) increase up to a local maximum and then goes asymptotically to zero. For second-order this local maximum occurs at value=1/2. and for first-order value=1/3. For second-order, below value=1/2. the up-wind finite difference approximation along the x direction, used to estimate the travel-time gradient, falls below 1/2 which may be physically incorrect. For first order, below value=1/3. the finite-difference approximation in the x direction falls below zero.

I still do not understand exactly what is happening. The solution may be to detect when the gradient at a point meets these conditions and abandon one of the directions of the stencil? I will be away from the computer until August 12th but can have a look when I get back.

[TomWinder]

Hi, have you made any further progress with this issue, or do you have any further advice how to handle it? I have recently updated to 2019.1.30 (from pip) from 0.0.9 and am now unable to use skfmm. Does this suggest the travel times I calculated with the previous version are likely to contain erroneous values? I'm happy to provide details of my specific use case if that would be useful. Thanks.

[jkfurtney]

Thanks for the note. If you could share some examples it would help to track down the problem. I have not looked at this recently but I will try to find some time soon. (I am a graduate of your department :) A student from the University of Exeter may be working on this module over the summer and they may be looking at this issue. @wmoebius

[jkfurtney]

OK, thanks to @f-fanni this is also now fixed in master.