Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

div by zero error in scikit-fmm extension velocity #15

Closed
conqueror5 opened this issue Jun 15, 2017 · 9 comments
Closed

div by zero error in scikit-fmm extension velocity #15

conqueror5 opened this issue Jun 15, 2017 · 9 comments
Labels

Comments

@conqueror5
Copy link

I am testing this package with random data and got this error message below:

terminate called after throwing an instance of 'std::runtime_error'
what(): div by zero error in scikit-fmm extension velocity

I am not sure whether there is any bugs or limitations of fast marching method. I attached my code.

Thanks
Code.zip

@jkfurtney
Copy link
Member

Thanks for getting in touch and for your testing effort. I have not seen this error before. It looks like a point which is exactly on the front is being finalized and the extension velocity cannot be calculated. I am not sure why this is happening, I will investigate in more detail. For now, this work around seems OK, could you test it?

In extension_velocity_marcher.cpp replace

  if (denominator != 0.0)
  {
    f_ext_[i] = numerator/denominator;
  }
  else
  {
    throw std::runtime_error(
          "div by zero error in scikit-fmm extension velocity");
  }

with

  if (denominator != 0.0)
  {
    f_ext_[i] = numerator/denominator;
  }
  else
  {
    double maxMagSpeed = 0;
    for (int dim=0; dim<dim_; dim++) {
      if (fabs(lspeed[dim]) > maxMagSpeed) {
        maxMagSpeed = lspeed[dim];
      }
    }
    f_ext_[i] = maxMagSpeed;
  }

The Adalsteinsson and Sethian paper has the details on this algorithm.

@conqueror5
Copy link
Author

It works for this example. Thank your for your help.

@jkfurtney
Copy link
Member

Here is the example that is not working:

data = {
	"phi": [
		[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
		[0,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,0],
		[0,-0.05,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.05,0],
		[0,-0.05,0,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-1.387778781e-17,-0.05,0],
		[0,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,0],
		[0,-0.05,0,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,4.163336342e-17,-0.05,0],
		[0,-0.05,0,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,0,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,2.636779683e-16,-0.05,0],
		[0,-0.05,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-6.938893904e-17,-0.05,0],
		[0,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,0],
		[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
	],
	"speed": [
		[11.45184936,11.88046603,13.93097355,15.93139249,15.8833956,14.84553532,13.6185673,12.42679836,11.31812304,10.29640191,9.358430712,8.497903081,7.708256266,6.983326657,6.317513336,5.70576199,5.1435114,4.626636758,4.151400612,3.714413028,3.312600057,2.943179047,2.603639393,2.291727521,2.005435111,1.742989795,1.502847699,1.283687324,1.084404378,0.9041071845,0.7421123969,0.5979407216,0.4713124028,0.3621422228,0.270533794,0.196772951,0.1413200935,0.1048014112,0.08799902395,0.09184022089,0.1173861727,0.1658207369,0.2384402707,0.3366457223,0.4619386639,0.6159233624,0.800317411,1.016973736,1.267916725,1.555393606,1.881938681,2.250426206,2.664064647,3.12570811,3.637418443,4.201174909,4.865898893,5.134755195,4.211308912,2.099396577,0.8102250364],
		[6.670346126,6.796495262,8.08806792,9.488340774,9.50478236,8.887207268,8.083054791,7.300475927,6.578794982,5.923831524,5.332620391,4.799544171,4.318858898,3.885246748,3.493941068,3.140703914,2.821772226,2.533802499,2.273822623,2.039191986,1.827568912,1.636884144,1.465319144,1.31128821,1.17342362,1.050563164,0.9417395949,0.8461716117,0.7632560729,0.6925611894,0.6338204799,0.586927288,0.5519296813,0.5290255578,0.5185578043,0.5210093776,0.5369982225,0.5672720077,0.6127027532,0.6742815571,0.7531137987,0.8504154095,0.9675110681,1.10583548,1.266939255,1.452501244,1.66434953,1.904493339,2.175167621,2.478889182,2.818516269,3.19726894,3.618634199,4.085579316,4.59957979,5.14635966,5.703471599,6.025970106,4.54286961,2.493210867,1.969980753],
		[1.82347227,1.735287094,1.629252084,1.630848852,1.61500077,1.479840409,1.279644864,1.087642323,0.9164628657,0.7692634732,0.644432858,0.5393962527,0.4516271379,0.3788446207,0.3190471898,0.2704946016,0.2316791724,0.2012979132,0.1782283676,0.161508263,0.1503184407,0.1439684509,0.1418842972,0.1435979124,0.1487380517,0.1570223712,0.1682505093,0.1822980443,0.1991112241,0.2187023874,0.2411460096,0.2665753097,0.2951793642,0.3272006759,0.3629331532,0.4027204653,0.4469547625,0.496075778,0.5505703676,0.6109726074,0.6778646413,0.7518785658,0.8336997562,0.9240721713,1.023806327,1.133790777,1.255008046,1.388555943,1.535674679,1.697778318,1.876484065,2.073610026,2.291092576,2.53057103,2.792709227,3.062913125,3.297050644,3.538802502,3.839664255,3.797256683,3.595070211],
		[1.400960471,1.267399181,0.6461617624,0.1079138162,0.0271989036,-0.009276508004,-0.0187279185,-0.02394268317,-0.02766901909,-0.03064419668,-0.03312135669,-0.0352148662,-0.03699144012,-0.03849796824,-0.03977126588,-0.04084191349,-0.04173600854,-0.04247609697,-0.04308174913,-0.04356996278,-0.04395546972,-0.04425098114,-0.04446738948,-0.04461393729,-0.0446983601,-0.04472700814,-0.04470495073,-0.04463606618,-0.04452311942,-0.04436782926,-0.04417092656,-0.04393220482,-0.04365056402,-0.04332404877,-0.04294988125,-0.04252448944,-0.04204353059,-0.04150190927,-0.04089378881,-0.0402125935,-0.03945099829,-0.03860090038,-0.03765336563,-0.03659853977,-0.03542551125,-0.0341221069,-0.03267459306,-0.0310672358,-0.02928163329,-0.02729563433,-0.02508140596,-0.02260143722,-0.01979895048,-0.01657063295,-0.01269204602,-0.00748676435,0.001822926996,0.05546429255,2.274156326,4.038789336,3.632126396],
		[0.6728483995,0.4966235082,0.1610875734,-0.02386456777,-0.04949632153,-0.04941923435,-0.04926075111,-0.04906091093,-0.04883053872,-0.04858101311,-0.04832127919,-0.04805778519,-0.04779521997,-0.0475371088,-0.04728617831,-0.04704456648,-0.04681395122,-0.0465956398,-0.04639063821,-0.04619970853,-0.04602341724,-0.04586217606,-0.04571627625,-0.04558591726,-0.04547123056,-0.04537229925,-0.04528917422,-0.04522188729,-0.0451704618,-0.04513492089,-0.04511529391,-0.04511162082,-0.04512395489,-0.04515236362,-0.04519692763,-0.04525773761,-0.04533488893,-0.04542847357,-0.04553856919,-0.04566522464,-0.04580844161,-0.04596815157,-0.04614418748,-0.04633624914,-0.04654386138,-0.04676632355,-0.04700264898,-0.04725149176,-0.04751105686,-0.0477789843,-0.04805218753,-0.04832660575,-0.04859680742,-0.04885540239,-0.0490923904,-0.04929559223,-0.04945218129,-0.04952383189,1.299192015,2.491535578,2.33515669],
		[0.2010596977,0.1506609227,0.02697122618,-0.04794332981,-0.04956536099,-0.04950032817,-0.04933721059,-0.04911892469,-0.04887255257,-0.04861289947,-0.04834779465,-0.04808181396,-0.0478182155,-0.04755966585,-0.04730844017,-0.0470664636,-0.04683532948,-0.04661632768,-0.04641048236,-0.046218593,-0.04604127366,-0.04587898825,-0.0457320811,-0.04560080286,-0.04548533228,-0.04538579432,-0.0453022751,-0.04523483413,-0.04518351412,-0.0451483487,-0.04512936815,-0.04512660333,-0.0451400879,-0.04516985873,-0.04521595451,-0.04527841249,-0.04535726308,-0.04545252224,-0.04556418134,-0.0456921943,-0.04583646155,-0.04599681059,-0.04617297265,-0.04636455496,-0.04657100814,-0.04679158849,-0.04702531529,-0.04727092445,-0.04752682078,-0.04779103046,-0.04806114427,-0.04833420278,-0.04860636933,-0.04887203155,-0.0491217619,-0.04933880616,-0.04949486607,-0.04955381603,0.5902321868,1.109538793,0.9890588414],
		[0.08096870911,0.07385511141,0.009852639914,-0.04830176362,-0.04956069614,-0.04948820999,-0.04932182255,-0.04910266879,-0.04885666254,-0.04859820007,-0.04833484851,-0.04807093696,-0.04780949903,-0.04755302097,-0.04730365087,-0.04706323769,-0.04683334187,-0.04661525404,-0.04641002379,-0.04621849276,-0.04604132749,-0.04587904949,-0.04573206185,-0.04560067203,-0.04548511132,-0.04538555116,-0.04530211674,-0.04523489816,-0.04518395937,-0.0451493452,-0.04513108643,-0.0451292032,-0.04514370672,-0.04517459927,-0.04522187256,-0.04528550435,-0.04536545323,-0.04546165165,-0.04557399685,-0.04570233985,-0.04584647224,-0.04600611075,-0.0461808794,-0.04637028923,-0.04657371533,-0.04679037142,-0.04701928246,-0.04725925649,-0.04750885759,-0.04776637996,-0.04802981093,-0.04829673072,-0.04856399655,-0.04882687078,-0.04907714707,-0.04930010701,-0.04947080869,-0.04955100671,0.2277990717,0.4853967311,0.4656325243],
		[0.1445421678,0.1368053297,0.06266536664,-0.02661682114,-0.04947557469,-0.04938176261,-0.04921536693,-0.04901335596,-0.04878438954,-0.04853870581,-0.04828439605,-0.04802713737,-0.0477709584,-0.04751887516,-0.04727326631,-0.04703606943,-0.04680888465,-0.04659303613,-0.04638961558,-0.04619951717,-0.04602346757,-0.04586205222,-0.04571573827,-0.04558489455,-0.04546980876,-0.04537070208,-0.0452877415,-0.04522104993,-0.04517071439,-0.04513679223,-0.04511931571,-0.04511829481,-0.04513371841,-0.04516555395,-0.04521374544,-0.04527820997,-0.04535883267,-0.04545546025,-0.04556789303,-0.04569587565,-0.04583908658,-0.04599712648,-0.04616950562,-0.04635563046,-0.04655478938,-0.04676613744,-0.04698867929,-0.04722124809,-0.04746247502,-0.04771073658,-0.04796405319,-0.0482198912,-0.04847480637,-0.04872393317,-0.04896061691,-0.0491780789,-0.04937221745,-0.0495084601,0.2176345573,0.633943082,0.7830633055],
		[0.3484777237,0.617460887,0.5291655426,0.1327816812,0.06101892167,0.01903577437,0.004379677977,-0.004771039515,-0.01161948951,-0.01713722218,-0.02170099158,-0.02551149492,-0.02870406725,-0.03138222762,-0.03362983769,-0.03551645301,-0.03710027024,-0.03843009729,-0.03954685358,-0.04048479402,-0.04127253819,-0.04193394398,-0.04248884876,-0.04295369511,-0.04334205462,-0.04366506159,-0.04393176682,-0.04414941991,-0.04432368717,-0.04445881036,-0.04455771018,-0.04462203671,-0.04465216742,-0.04464715165,-0.04460459879,-0.04452050543,-0.04438901481,-0.04420209972,-0.04394915791,-0.04361650676,-0.04318676227,-0.04263808549,-0.04194327888,-0.04106871505,-0.0399730815,-0.03860592727,-0.03690599842,-0.03479934508,-0.03219715788,-0.02899320414,-0.02506047553,-0.02024574392,-0.01435830143,-0.00713629306,0.001858536507,0.01375986518,0.03143151673,0.08135764101,0.8104829951,1.455284408,1.410104882],
		[0.3549682693,0.9170393495,2.391263054,3.777964661,4.157454248,3.787378189,3.253506372,2.76392184,2.334022702,1.965084181,1.65121391,1.384959523,1.159268353,0.9679267404,0.8056049537,0.6677835091,0.5506526884,0.4510129847,0.3661848745,0.2939296571,0.2323810492,0.1799866158,0.1354580121,0.09772903604,0.06592057877,0.03931165781,0.01731582759,-0.0005376219893,-0.01461816196,-0.02520458882,-0.0324919159,-0.03659474351,-0.03754714759,-0.03529901078,-0.02970859424,-0.02053102238,-0.00740221003,0.01018238549,0.03289495894,0.06160942388,0.09744342015,0.1418094509,0.1964764131,0.263642793,0.346022703,0.4469457038,0.570470917,0.7215151629,0.9059934551,1.13096682,1.404786064,1.737190941,2.139314765,2.622946803,3.199441285,3.868694272,4.692876853,5.123769267,4.11750588,2.460239563,1.772368201],
		[0.5122036578,1.111692571,4.814691278,8.864074419,9.455304413,8.497019134,7.317599258,6.232429749,5.283653852,4.469041582,3.774465976,3.183604169,2.681224665,2.253937136,1.890249492,1.580402923,1.316157682,1.090582174,0.8978616632,0.7331300774,0.5923243322,0.4720593425,0.3695216121,0.2823793553,0.208707255,0.1469241733,0.09574235822,0.05412691844,0.02126456652,-0.003459136088,-0.02047457322,-0.03003521563,-0.03221718788,-0.02691109865,-0.01380571165,0.007637231791,0.03821802042,0.07904509899,0.1315939205,0.1977804509,0.2800515477,0.3814946799,0.5059696076,0.6582646454,0.8442799219,1.071239538,1.347933577,1.684989249,2.095167501,2.593674496,3.198465296,3.93045705,4.813570673,5.873442898,7.136760586,8.603921568,10.33486172,11.6700284,7.97361342,2.821032776,1.750295625],
		[0.8981998045,1.510741814,3.457929792,5.221793674,5.358778003,4.728586308,4.068364103,3.463602339,2.93784454,2.48661732,2.101311656,1.77285755,1.492941385,1.25428324,1.05063706,0.8766940953,0.7279662353,0.6006726844,0.4916373914,0.3981987927,0.3181315149,0.2495790832,0.1909965326,0.1411018526,0.09883527464,0.06332551762,0.0338622321,0.009873995351,-0.009088665591,-0.02336562246,-0.03319515153,-0.03871775534,-0.03997583889,-0.03690915367,-0.02934578557,-0.01698832078,0.0006053282487,0.02404711623,0.05415092502,0.09197250417,0.1388584597,0.196505587,0.2670319188,0.3530608597,0.4578196593,0.5852531976,0.7401535432,0.9283048542,1.156641652,1.433415051,1.768355996,2.172794652,2.659708554,3.243205622,3.939054327,4.748887929,5.673335356,6.627542937,5.906527178,3.083222187,1.442373667],
		[0.6075922344,0.8577360822,0.5548156836,-0.023943917,-0.04960909597,-0.04952125463,-0.04939727743,-0.04925554672,-0.04909696628,-0.04892528563,-0.04874521655,-0.04856098233,-0.04837594335,-0.04819270105,-0.04801328902,-0.04783933306,-0.04767216075,-0.04751287148,-0.04736238157,-0.04722145491,-0.04709072559,-0.04697071606,-0.04686185241,-0.04676447759,-0.04667886301,-0.04660521851,-0.04654370099,-0.04649442156,-0.0464574514,-0.04643282634,-0.04642055024,-0.04642059711,-0.04643291224,-0.04645741224,-0.04649398411,-0.04654248348,-0.04660273215,-0.04667451497,-0.04675757639,-0.04685161675,-0.04695628864,-0.0470711935,-0.04719587867,-0.04732983489,-0.04747249396,-0.04762322552,-0.0477813307,-0.04794602814,-0.04811642484,-0.04829146041,-0.04846981197,-0.04864975667,-0.04882903109,-0.04900484032,-0.04917444237,-0.04933639824,-0.04949074194,-0.04961174316,1.730371126,2.304935266,1.099470864],
		[0.08541049203,0.04212572542,-0.02456997689,-0.04881200858,-0.0496309446,-0.04957118772,-0.04945648119,-0.04930884225,-0.04913944176,-0.04895723641,-0.04876868203,-0.04857805988,-0.0483882816,-0.04820148695,-0.04801937301,-0.04784334956,-0.04767460864,-0.04751415838,-0.04736284478,-0.04722136999,-0.04709030983,-0.04697013063,-0.04686120499,-0.04676382622,-0.04667822105,-0.04660456074,-0.04654297052,-0.04649353746,-0.04645631686,-0.0464313372,-0.04641860386,-0.04641810154,-0.04642979548,-0.04645363167,-0.04648953581,-0.04653741143,-0.04659713689,-0.04666856159,-0.04675150128,-0.04684573272,-0.0469509877,-0.04706694671,-0.04719323227,-0.04732940231,-0.04747494373,-0.04762926614,-0.04779169546,-0.04796146558,-0.04813770244,-0.04831938762,-0.0485052722,-0.04869368674,-0.04888216441,-0.04906681795,-0.04924158223,-0.04939794263,-0.04952291842,-0.04958703342,0.4709075762,0.8920488717,0.7926840498],
		[-0.03315760667,-0.02764392508,-0.03541226377,-0.0491583532,-0.04962095017,-0.04958135369,-0.04947687328,-0.04932753206,-0.04915277637,-0.04896515396,-0.04877207962,-0.04857789572,-0.04838534294,-0.04819638751,-0.04801260283,-0.04783531848,-0.04766567102,-0.04750462283,-0.04735297527,-0.04721138368,-0.04708037463,-0.04696036418,-0.04685167548,-0.04675455499,-0.04666918694,-0.04659570569,-0.04653420628,-0.04648475305,-0.04644738671,-0.04642212978,-0.04640899064,-0.04640796624,-0.0464190435,-0.04644219951,-0.04647740043,-0.04652459924,-0.04658373217,-0.04665471393,-0.0467374316,-0.04683173731,-0.04693743947,-0.04705429289,-0.04718198765,-0.04732013698,-0.04746826464,-0.04762579221,-0.04779202691,-0.04796614945,-0.04814719806,-0.04833403498,-0.04852526015,-0.04871899485,-0.04891239361,-0.04910068353,-0.04927558475,-0.04942322139,-0.04952291196,-0.04955364601,0.2846575114,0.8707168826,1.122569206],
		[-0.00188915534,-0.003813862766,-0.02598423053,-0.04791583193,-0.04960045167,-0.0495564095,-0.04944698334,-0.04929896024,-0.04912760688,-0.04894287627,-0.0487515108,-0.04855800972,-0.04836549666,-0.0481762881,-0.04799220307,-0.04781471628,-0.04764503426,-0.04748413903,-0.04733282008,-0.04719170216,-0.04706127114,-0.04694189776,-0.04683385917,-0.04673735782,-0.04665253777,-0.04657949855,-0.04651830671,-0.04646900532,-0.0464316216,-0.04640617289,-0.046392671,-0.04639112523,-0.04640154395,-0.04642393477,-0.04645830341,-0.04650465108,-0.04656297039,-0.04663323956,-0.04671541503,-0.04680942197,-0.04691514291,-0.04703240407,-0.04716095927,-0.04730047138,-0.04745049122,-0.04761043367,-0.04777955094,-0.04795690146,-0.04814131013,-0.04833130875,-0.04852503063,-0.04872000962,-0.04891280369,-0.04909834988,-0.04926903443,-0.04941345768,-0.04951373679,-0.04954438288,0.5482434556,1.896183799,2.646341997],
		[0.2002175433,0.1912373017,0.09196593107,-0.02392854761,-0.04953133654,-0.04948239159,-0.04937065577,-0.04923161984,-0.04907087396,-0.04889463802,-0.04870900464,-0.04851902535,-0.04832861372,-0.04814073942,-0.04795765825,-0.04778109764,-0.04761238902,-0.04745255983,-0.04730239892,-0.04716250551,-0.04703332749,-0.04691519267,-0.04680833461,-0.04671291429,-0.04662903821,-0.04655677347,-0.0464961602,-0.04644722171,-0.04640997263,-0.04638442524,-0.04637059421,-0.04636849982,-0.04637816969,-0.04639963913,-0.04643294995,-0.04647814764,-0.04653527695,-0.04660437537,-0.04668546457,-0.04677853939,-0.04688355406,-0.04700040518,-0.04712891122,-0.04726878758,-0.04741961665,-0.04758081124,-0.04775156942,-0.04793081676,-0.04811712997,-0.04830863211,-0.04850284784,-0.04869651171,-0.04888535237,-0.04906396532,-0.0492260365,-0.04936549672,-0.04947507732,-0.04952304225,1.303274937,3.936989477,5.217904262],
		[0.4486086894,0.640762828,0.4920788008,0.1028941073,0.02583280023,-0.008917565278,-0.01799359874,-0.0231292393,-0.0268951016,-0.02997236032,-0.03258519926,-0.0348313785,-0.03676792109,-0.03843573003,-0.03986760154,-0.04109122644,-0.04213056256,-0.04300658344,-0.04373774528,-0.04434030699,-0.04482856595,-0.04521504178,-0.0455106259,-0.0457247073,-0.04586528108,-0.04593904366,-0.04595147766,-0.04590692824,-0.0458086725,-0.04565898289,-0.04545918576,-0.04520971546,-0.04491016489,-0.04455933271,-0.04415526764,-0.0436953097,-0.04317612824,-0.04259375571,-0.0419436159,-0.04122054406,-0.04041879575,-0.03953203936,-0.03855332579,-0.03747502644,-0.03628872697,-0.03498505982,-0.03355344982,-0.03198173191,-0.03025556577,-0.02835747982,-0.02626509298,-0.02394703634,-0.02135163483,-0.01837042606,-0.01472268432,-0.009448051463,0.00180676815,0.08140788664,2.372375428,6.102216831,7.6126698],
		[0.8122135474,1.063671762,1.442337212,1.623889942,1.654390605,1.532924648,1.335323758,1.142528355,0.9679423557,0.8155049786,0.6843426256,0.5724507762,0.4777191966,0.3981642894,0.3319939606,0.2776124353,0.2336078335,0.1987359458,0.1719046751,0.1521602755,0.1386754469,0.130739032,0.1277470122,0.1291945416,0.1346688139,0.1438426129,0.1564684373,0.1723731231,0.1914529092,0.2136689022,0.2390429067,0.2676535901,0.299632953,0.3351630828,0.3744731697,0.4178367813,0.465569404,0.5180262903,0.5756006884,0.6387225741,0.7078580724,0.7835098277,0.8662186668,0.9565669898,1.055184413,1.162756236,1.280035266,1.407857222,1.547159024,1.698996518,1.864551491,2.045090636,2.241806747,2.455208044,2.683782719,2.908398487,3.072360103,3.206218387,3.641257682,5.743255038,7.487012806],
		[4.799723978,5.668785456,7.759713718,9.23798731,9.303146352,8.733942039,7.967617099,7.213592242,6.511925303,5.870184571,5.28713913,4.758606071,4.279951764,3.846703161,3.454734977,3.10029614,2.779991751,2.49075381,2.229812556,1.994671571,1.783087046,1.593050661,1.42277537,1.270683416,1.135396034,1.015724416,0.9106616142,0.8193751325,0.7412000234,0.6756323306,0.6223227499,0.5810703914,0.5518165359,0.5346382924,0.5297420728,0.5374568325,0.5582270602,0.5926055649,0.6412461893,0.7048966878,0.7843921499,0.8806495167,0.994663941,1.127507954,1.280334636,1.454386148,1.651009064,1.871677558,2.11802426,2.391874478,2.695269744,3.030425302,3.399510937,3.803574754,4.24068802,4.689174361,5.099799182,5.205761261,3.924302031,3.765180529,4.805009321],
		[8.502338763,10.1336972,13.42210996,15.43502472,15.39865245,14.43362972,13.27792196,12.1450601,11.08321226,10.09827969,9.189105504,8.351142862,7.579235281,6.86832537,6.213683631,5.61094594,5.056099736,4.54545654,4.075624896,3.643487528,3.246183253,2.881093013,2.545829182,2.238227301,1.95633955,1.698429406,1.462967036,1.2486251,1.05427468,0.8789811152,0.7219995508,0.5827700286,0.4609119596,0.3562178382,0.2686460688,0.1983128133,0.1454828072,0.1105591639,0.09407227802,0.09666806527,0.1190959325,0.1621970683,0.2268938743,0.3141816146,0.4251236442,0.5608518411,0.7225740654,0.9115904231,1.12931948,1.377333029,1.657392548,1.971454199,2.321570578,2.709032664,3.133562016,3.591986874,4.107934478,4.211401594,3.497573519,3.036945462,3.114693481]
	],
	"dx": 0.05
}
import skfmm
skfmm.extension_velocities(data["phi"], data["speed"], dx=data["dx"])

@jkfurtney
Copy link
Member

Interestingly, if you omit the dx keyword argument and let it default to 1.0 there is no crash??

@jkfurtney
Copy link
Member

jkfurtney commented Feb 4, 2019

Here is a simplified case that shows the problem:

import numpy as np
import skfmm

phi = np.array([[  0.00000000e+00],
       [  5.00000000e-02],
       [  1.00000000e-01],
       [  1.50000000e-01],
       [  1.50000000e-01],
       [  1.00000000e-01],
       [  5.00000000e-02],
       [ -1.38777878e-17]])

speed = np.array([[ 0.2704946 ],
       [-0.04084191],
       [-0.04704457],
       [-0.04706646],
       [-0.04706324],
       [-0.04703607],
       [-0.03551645],
       [ 0.66778351]])
dx = 0.05
d, f_ext = skfmm.extension_velocities(phi, speed, dx=dx)
  • If dx = 1.0, the crash does not happen
  • phi and speed here are 2D arrays but only have one column of data. If you run it as a 1D problem it works
  • Possibly the negative speed values are a problem?

murraycutforth added a commit to murraycutforth/scikit-fmm that referenced this issue Sep 21, 2020
@murraycutforth
Copy link
Contributor

Hi, I've also encountered a div by zero error in extension_velocities. I'm pretty sure that it's the same problem as described here. I looked into the exact values of phi in surrounding cells at the point when the error occurred in my case, and found that it reduced to the following 1D problem as the neighbouring cells in the other two dimensions were not frozen:

image

The issue is happening when the wave front propagating outward arrives at a point at the same time from both directions, so d, and ldistance[dim] are set to zero.

for (int dim=0; dim<dim_; dim++)
{
lspeed[dim]=0;
ldistance[dim]=0;
for (int j=-1; j<2; j+=2) // each direction
{
int naddr = _getN(i,dim,j,Mask);
if (naddr!=-1 && flag_[naddr]==Frozen)
{
//determine which direction, in this dimension, is nearest to
//the front. Calculate the distance to front in this direction
double d = distance_[i] - distance_[naddr];
if (ldistance[dim]==0 || ldistance[dim]>d)
{
ldistance[dim] = d;
lspeed[dim] = f_ext_[naddr];
}
}
} // for each direction
} // for each dimension

I think a viable fix is just to check that d is not zero before updating ldistance and lspeed, like this:

https://github.com/murraycutforth/scikit-fmm/blob/b1cf74c5ed14a30c3b1432c9cec38399ebb531ff/skfmm/extension_velocity_marcher.cpp#L122-L126

This runs successfully on my case, as well as on the simplified test case posted above. The extension velocity example script also remains correct. This fix means that the velocity is extended from the neighbouring cell which phi was calculated from, so the extension velocity and phi are consistent with each other (e.g. in my hand-drawn example above, v is extended from the cell on the right hand side, rather than the left hand side).

Do you think this is the right approach @jkfurtney ?

@jkfurtney
Copy link
Member

Thanks for your investigation here and for your intuitive explanation. What you describe makes sense, could you put this in a PR? I will merge it and give it a few more tests.

@jkfurtney
Copy link
Member

Is that fountain pen ink?

@murraycutforth
Copy link
Contributor

Sure, just sent the PR. I don't think it should change the result of the function since any d > eps is treated the same as before, but would be great if you can run some regression tests to double check.

Haha yeah almost, I used to use a fountain pen but switched to one of those uni-ball pens now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants