assertion issue in Smooth.h

[cmaurini]

I got a error for the assertion at line 838 in Smooth.h: assert(mag!=0); if I comment this line the code seems to work correctly, but then there should a div by 0 in line 841. Could you check this? I got this on macosx and clang.

[croci]

Would you be able to provide a minimum working example that reproduces the error?

Then I am out, but maybe @taupalosaurus can help

[cmaurini]

With this code

from dolfin import *
from adaptivity import *
mesh = UnitSquareMesh(5,5)
M  = mesh_metric(mesh)
Mp = refine_metric(M, 2)
new_mesh, boundary_tags = adapt(Mp)

I get

target mesh has 87 nodes
Colinearity detection ...
Beginning PRAgMaTIc adapt
Initialising PRAgMaTIc ...
Setting surface ...
Setting metric tensor field ...
Entering adapt ...
Assertion failed: (mag!=0), function optimisation_linf_2d_kernel, file /Users/maurini/Documents/MyDoc/WIP/codes/pragmatic/include/Smooth.h, line 838.
[CMs-MacBook-Pro:99084] *** Process received signal ***
[CMs-MacBook-Pro:99084] Signal: Abort trap: 6 (6)
[CMs-MacBook-Pro:99084] Signal code:  (0)
[CMs-MacBook-Pro:99084] [ 0] 0   libsystem_platform.dylib            0x00007fff54ad7f5a _sigtramp + 26
[CMs-MacBook-Pro:99084] [ 1] 0   ???                                 0x0000000000000020 0x0 + 32
[CMs-MacBook-Pro:99084] [ 2] 0   libsystem_c.dylib                   0x00007fff5490330a abort + 127
[CMs-MacBook-Pro:99084] [ 3] 0   libsystem_c.dylib                   0x00007fff548cb360 basename_r + 0
[CMs-MacBook-Pro:99084] [ 4] 0   libpragmatic.dylib                  0x000000011b479cae _ZN6SmoothIdLi2EE27optimisation_linf_2d_kernelEi + 2526
[CMs-MacBook-Pro:99084] [ 5] 0   libpragmatic.dylib                  0x000000011b4792bb _ZN6SmoothIdLi2EE24optimisation_linf_kernelEi + 27
[CMs-MacBook-Pro:99084] [ 6] 0   libpragmatic.dylib                  0x000000011b3717cb _ZN6SmoothIdLi2EE17optimisation_linfEid + 5115
[CMs-MacBook-Pro:99084] [ 7] 0   libpragmatic.dylib                  0x000000011b35768a pragmatic_adapt + 874
[CMs-MacBook-Pro:99084] [ 8] 0   _ctypes.cpython-36m-darwin.so       0x000000010f50849f ffi_call_unix64 + 79
[CMs-MacBook-Pro:99084] [ 9] 0   ???                                 0x00007ffee15c5c50 0x0 + 140732679347280
[CMs-MacBook-Pro:99084] *** End of error message ***
Abort trap: 6

Maybe related to uniform metric? The python tests pass otherwise.

[taupalosaurus]

Hi @cmaurini, I am not sure what is going on there. I am going to try to reproduce the case. Are you using the latest version of pragmatic ? Thanks

[croci]

Ok, the above script runs on my local FEniCS installation on my laptop and in a FEniCS/stable and FEniCS/dev docker container in both python2.7 and python3. So it could be an issue with either macosx, clang or your installation in general. Sorry for not being able to help more.

[taupalosaurus]

@cmaurini were you using the latest version of pragmatic ? I fixed a bug in the smoothing code in september... Otherwise, can you please output from FEniCS the mesh and metric that make the code fail ? I'd just like to check we aren't missing the obvious before bringing out the big guns! Thanks

[cmaurini]

Yes, I use last version.

thanks @taupalosaurus for looking at this.

To me the explication is simple: the metric we give is uniform, i.e. a tensor constant in space, so I suppose that the grad_w given by lipnikov_grad(loc, x0, x1, x2, x3, m0, grad_w); is null and then is norm double mag = sqrt(grad_w[0]*grad_w[0] + grad_w[1]*grad_w[1] + grad_w[2]*grad_w[2]); is null. So assert(std::isnormal(mag)); fails. As it should (http://en.cppreference.com/w/cpp/numeric/math/isnormal).

This can be detected or not as of function of the compiler of the options I suppose, but it is a bug in the code to me: it does not accept uniform metrics as input. But I admit that I am completely ignorant on the detail of the code there.

I really I am not an expert of the used algorithm, but I suppose that a easy fix could be to define mag in a regularized way, as 'double mag = sqrt(grad_w[0]grad_w[0] + grad_w[1]grad_w[1] + grad_w[2]*grad_w[2]+10e-14);' but maybe it can give some numerical issue. Otherwise one could do something different than search[i] = grad_w[i]/mag; for uniform meshes?

[taupalosaurus]

The algorithm is more complicated: what lipnikov_grad returns is the variation of the quality (wrt to metric m0) of the element x0, x1, x2 when you move x0. The fact that the metric is uniform doesn't change the fact this quantity should not be null. Uniform metrics are heavily tested, as it is the simplest test possible, and given that Matteo cannot reproduce the bug, I fear that it is much more vicious... Which is why I wanted to check your input metric first ( @croci fixed something this morning, could it be related ?) Thanks,

[croci]

What I fixed this morning is unrelated and I ran all the tests as well as the basic example of the issue after the change both in python2 and 3 and it runs.

[cmaurini]

@taupalosaurus This should give you the information you need. I printed out in a pdf containing images, codes, and data about the metric. The mag gives numbers of the order of 10^-8, so I am not surprised by the assertion error. Maybe what happens running it depends on the tolerance of this assertion test on different compilers if it gives an error or not, but it is almost zero in any case.

test.pdf

[taupalosaurus]

thanks @cmaurini The good news is I managed to install fenics on my macbook and I can try to understand what is going on now!

[taupalosaurus]

So I was able to reproduce this bug, in this case mag is indeed 0 (to the 15th digit). What that means is we are in a very rare configuration where numerical precision results in the lipnikov functional of the slightly perturbed triangles are the same within machine precision. I don't see anything wrong in this part of the code, and I don't really see why this is catched by an assert. In special configurations where the metric is the one of the rectangle triangle, it is likely to happen again. @ggorman Do you remember why you put an assert here, instead of catching the 0 and just doing nothing ?

[taupalosaurus]

I have merged a fix: if the gradient is null, it means the vertex is already in optimal position so I don't try to smooth it. Hopefully that is enough. Thanks @cmaurini for the issue, please do let me know if there is something else wrong!