Point interpolation is (too?) slow

[stephankramer]

Despite @ReubenHill 's improvements in 81644ff I have found point interpolation for large point sets to still be fairly slow. I in particular noticed this comparing mesh-to-mesh interpolation in an adaptive mesh simulation with the same setup in fluidity. In fluidity the time spent on interpolation is pretty much negligible whereas in my firedrake setup it dominates by an order of magnitude - taking the time of a dozen Stokes timesteps, making adaptivity pointless. I know @jrmaddison and @pefarrell have in the past spent a lot of time optimising the node/intersection finder in fluidity and much of this has been transferred to libsupermesh. So this brought me to the idea of just trying to use the rtree node locator used internally by libsupermesh and indeed doing so it seems I can bring down the cost of interpolation by a factor of 20.

Take the following test example:

from firedrake import *

nx, ny = 100, 100
mesh = UnitSquareMesh(nx, ny, quadrilateral=False)
Q = FunctionSpace(mesh, "CG", 2)
u1 = Function(Q)
u2 = Function(Q)

def interpolate_between(source, target):
    from firedrake.mg.utils import physical_node_locations
    target_fs = target.function_space()
    target_locations = physical_node_locations(target_fs).dat.data[:]
    vom = VertexOnlyMesh(source.function_space().mesh(), target_locations, redundant=False)
    if isinstance(target_fs.ufl_element(), VectorElement):
        vom_space = VectorFunctionSpace(vom, "DG", 0)
    else:
        vom_space = FunctionSpace(vom, "DG", 0)
    target_vom_f = interpolate(source, vom_space)
    target.dat.data[:] = target_vom_f.dat.data

interpolate_between(u1, u2)

This interpolates to ~ 40,000 nodes and takes about 2 minutes on master whereas using the rtree configuration of libsupermesh it takes less than 5 seconds. To reproduce the latter you need https://bitbucket.org/stephankramer/libsupermesh/branch/skramer/expose-rtree-c-interface to expose the c interface to libspatialindex's rtree as it's internally used in libsupermesh and use the skramer/libsupermesh-rtree branch of Firedrake.

To be honest, I don't really understand why firedrake's use of libspatialindex's rtree is so much slower. Index_Intersects_id just creates a Region and calls the intersectsWithQuery() method on the rtree index for every queried point, which is pretty much what is done in libsupermesh_query_rtree() except that it creates a Point. I tried writing my own version of Index_Intersects_id that still uses firedrake's created rtree object but uses a Point instead, but that didn't seem to make any difference. So it must be something in the configuration of the rtree index that is different.

[stephankramer]

I have profiled the above test case using py-spy, and this is what I get on master:

and this is what I get using libsupermesh's rtree config:

[jrmaddison]

Are the costs due to the interpolation, or building the tree? libsupermesh uses bulk loading, and afaict Firedrake is adding nodes one at a time using Index_InsertData. libsupermesh also avoids the use of disk caching by libspatialindex.

[stephankramer]

Just querying the tree: it spends all its time in Index_Intersects_id (90%) and just 1.24% in Index_InsertData. That disk caching sounds interesting: is that just a configuration option from libspatialindex or requires some additional code?

[jrmaddison]

I've only seen disk caching when bulk loading.

[ReubenHill]

Some thoughts:

• Do we have the same bounding boxes in each case? And do we follow the same code paths (e.g. using locate_cell) in both cases?
• Do tests all pass with libsupermesh's rtree? If the only difference is using the one in libsupermesh then maybe we can swap it in?

[stephankramer]

Ah, I figured out I've been completely barking up the wrong tree (!) :
My issue is actually caused by the change in default tolerance, the tolerance that the boundary boxes of all individual elements are increased by. The default is now 1.0 which is not a very good default for a unit square :-) as it makes all bounding boxes the size of the domain, and so all points will return all elements as possible candidates. If I set the tolerance correctly for my case to something small I get roughly the same performance with firedrake's rtree.

I guess it would be a good idea to give UnitSquareMesh a better default for this, but more generally I've never quite liked having to use an absolute tolerance for this. I don't think this works very well for multi-scale meshes: if I have cells that are 1 km and cells that are 5m what tolerance should I use? Would it not be much better to simply expand the bounding boxes with a relative tolerance, i.e. relative to the extents of the individual elements? Fluidity simply added 10% - if your mesh gradates reasonably well, that doesn't really add to the number of possible candidates much, and should be plenty to keep points on the boundary inside. Or perhaps have both an absolute and relative tolerance, but setting the default absolute tolerance to 0.

[ReubenHill]

The new tolerance is cell-size relative so 1.0 does make sense. So if you have a 100x100 square mesh of quads the tolerance in global space is 0.01. EDIT: for triangles, as used here, it's the same order of magnitude

[dham]

Ah, but that's currently bugged pending Kaito's fix. This is the bit where you added a relative tolerance to an absolute value. Current master effectively uses an absolute tolerance of 1.0. We should get Kaito's thing merged.

[ReubenHill]

Kaito's fix extends the bounding boxes to another dimension when the cells are immersed manifolds. The diff (which only does this when there is a tolerance which I think he needs to fix) is

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index fc3aec077..1eb5d5b82 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2014,8 +2014,10 @@ values from f.)"""
         # tolerance. Note that if tolerance is too small it might not actually
         # change the value!
         if hasattr(self, "tolerance") and self.tolerance is not None:
-            coords_min -= self.tolerance
-            coords_max += self.tolerance
+            coords_mid = (coords_max + coords_min)/2
+            d = np.max(coords_max - coords_min, axis=1)[:, None]
+            coords_min = coords_mid - (self.tolerance + 0.5)*d
+            coords_max = coords_mid + (self.tolerance + 0.5)*d
 
         # Build spatial index
         return spatialindex.from_regions(coords_min, coords_max)

Can you try this out and see if it fixes the problem with the original tolerance?

[ReubenHill]

Or indeed this slightly simpler fix which is along the lines of what we discussed on slack just now

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index fc3aec077..aea3c210f 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2014,8 +2014,9 @@ values from f.)"""
         # tolerance. Note that if tolerance is too small it might not actually
         # change the value!
         if hasattr(self, "tolerance") and self.tolerance is not None:
-            coords_min -= self.tolerance
-            coords_max += self.tolerance
+            coords_diff = (coords_max - coords_min)
+            coords_min -= self.tolerance*coords_diff
+            coords_max += self.tolerance*coords_diff

         # Build spatial index
         return spatialindex.from_regions(coords_min, coords_max)

[ReubenHill]

Opened #2909 which fixes this. If it passes tests then it should be merged.

[stephankramer]

Thanks @ReubenHill that fixes the performance bug for me. The example now runs in 6s rather than 2 min

[stephankramer]

Fixed by #2909