Broken locate_cells_ref_coords_and_dists?

[rubyturner]

I am having some extremely bizarre problems with the locate_cells_ref_coords_and_dists mesh method which I noticed when trying to use VertexOnlyMesh. It seems to be specific to the installation on my Mac - I haven't been able to reproduce the problem anywhere else.

For example, something like this:

from firedrake import *
import numpy as np

mesh = UnitSquareMesh(3,3, quadrilateral=True)
coords = np.array([[0,0],[0.9,0.1],[0.5,0.5],[0.2, 0.8]])
print(mesh.locate_cells_ref_coords_and_dists(coords))

Will sometimes give:

(array([          0,  1073741824,           0, -1074790400], dtype=int32), 
 array([[0.e+000, 0.e+000],
       [0.e+000, 5.e-324],
       [5.e-324, 0.e+000],
       [5.e-324, 5.e-324]]), 
 array([0.        , 0.33333333, 0.66666667, 1.        ]))

or otherwise,

(array([3, 0, 7, 0], dtype=int32), 
 array([[0.e+000, 0.e+000],
       [0.e+000, 5.e-324],
       [5.e-324, 0.e+000],
       [5.e-324, 5.e-324]]), 
 array([0.        , 0.33333333, 0.66666667, 1.        ]))

rather than the expected:

(array([0, 3, 4, 5], dtype=int32),
 array([[0. , 0. ],
        [0.7, 0.3],
        [0.5, 0.5],
        [0.6, 0.4]]),
 array([0., 0., 0., 0.]))

I can't find any other problems (only the vertexonly tests fail), and I cannot pinpoint the issue any further.

I have tried updating and re-installing firedrake to no avail. And there are no warnings or issues with the install at all.

Has anyone else run in to this problem before? I would appreciate any pointers as to how to resolve or any next steps.

[dham]

I can reproduce this:

In [1]: from firedrake import *
   ...: import numpy as np
   ...: 
   ...: mesh = UnitSquareMesh(3,3, quadrilateral=True)
   ...: coords = coords = np.array([[0,0],[0.9,0.1],[0.5,0.5],[0.2, 0.8]])
   ...: print(mesh.locate_cells_ref_coords_and_dists(coords))
(array([3, 0, 7, 0], dtype=int32), array([[0., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.]]), array([0.        , 0.33333333, 0.66666667, 1.        ]))

I also get a bunch of test fails including a load of out of bounds errors. Looks like we have a memory access bug on Mac.

[dham]

Yeah. No beans.

What's happening is that the reference cells, reference cell coordinates and distances from cells are not getting updated by the call to self._c_locator in mesh.locate_cells_ref_coords_and_dists. The values above are just whatever starts out in those memory locations when numpy.empty creates the arrays. This explains the seemingly stochastic behaviour that Ruby observed.

My current hypothesis is that somehow the pointer that the c function is getting is not the right one.

[wence-]

The ctypes "pointer types" in that function are hard-coded as double and int. It is possible (though unlikely) that those are not matched up with the actual datatypes being passed in? Ah, no. There's a mismatch in the way the pointer is called and the arguments it is declared with:

Call-site:

firedrake/firedrake/mesh.py

Lines 2142 to 2147 in 7eb993f

	self._c_locator(tolerance=tolerance)(self.coordinates._ctypes,
	xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
	Xs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
	ref_cell_dists_l1.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
	cells.ctypes.data_as(ctypes.POINTER(ctypes.c_int)),
	npoints)

Which matches the C function signature:

firedrake/firedrake/mesh.py

Line 2162 in 7eb993f

int locator(struct Function *f, double *x, double *X, double *ref_cell_dists_l1, int *cells, size_t npoints)

But, ctypes-declaration site:

firedrake/firedrake/mesh.py

Lines 2191 to 2194 in 7eb993f

	locator.argtypes = [ctypes.POINTER(function._CFunction),
	ctypes.POINTER(ctypes.c_double),
	ctypes.POINTER(ctypes.c_double),
	ctypes.POINTER(ctypes.c_double)]

So you're missing two arguments there.

[wence-]

Aside, this code should validate that the passed in arrays and coordinate arrays have dtypes that match the C function (otherwise you'll also get garbage); whatever.ctypes.data_as(...) does not do any numpy astype casting, it just says "interpret this array as a pointer to foo".

[dham]

I claim #3003 fixes this. Thanks for the very clear bug report @rubyturner.

[dham]

The fix has now landed in master so hopefully running firedrake-update will solve everything.

[ReubenHill]

Weird that this slipped through given that I run my vertex only mesh tests on macos and our CI passes. How is it possible that this wasn't noticed?

Aside, this code should validate that the passed in arrays and coordinate arrays have dtypes that match the C function (otherwise you'll also get garbage); whatever.ctypes.data_as(...) does not do any numpy astype casting, it just says "interpret this array as a pointer to foo".

I wonder where the best location for this would be... Is this pattern only used here or are there other places we should be checking?

[rubyturner]

Thanks for the fix! I can confirm it is all working now for me. Does anyone have any insight as to why this bug only caused issues on some machines? It was working fine on a different mac.

[wence-]

I wonder where the best location for this would be... Is this pattern only used here or are there other places we should be checking?

Well, any time you use ctypes to provide a signature for a dlsym'ed symbol you need to ensure that: (1) that signature is correct; (2) the call of that symbol is correct.

[dham]

Thanks for the fix! I can confirm it is all working now for me. Does anyone have any insight as to why this bug only caused issues on some machines? It was working fine on a different mac.

This is an error in how we were calling a C function. C is a very low-level language, and this sort of error results in accessing memory at the wrong locations. Unlike Python, where analogous errors will immediately raise an exception, in C the effect of this sort of bug is sensitively dependent on the precise layout of the memory in the program concerned. This depends on hardware, operating system, compiler version, library versions.... There are any number of seemingly innocuous changes that can change the behaviour of a bug like this. Indeed, the question that presents itself when we find an error of this sort isn't really "why did that fail now?", it's "how on Earth did we get away with this for so long!"