FEniCS / FEniCS/dolfinx

Expose template type of `compute_distance_gjk` to Python users

Open
#3,836 1 comment 2 reactions 0 assignees View on GitHub
enhancement performance
Dominant language
C++
Stars
1.2k
Forks
261
Avg merge
1d 15h
Merged PRs (30d)
65

Description

For spatial searches, the function `compute_distance_gjk` is templated over scalar types and defaults to using `boost::multiprecision::cpp_bin_float_double_extended>`:

https://github.com/FEniCS/dolfinx/blob/3817e5174080d233aebdb7e87793ba7696548616/cpp/dolfinx/geometry/gjk.h#L306-L309

This function is heavily relied upon by `determine_point_ownership` via `compute_first_colliding_cell`. The call inside of `compute_first_colliding_cell` looks like this:

https://github.com/FEniCS/dolfinx/blob/3817e5174080d233aebdb7e87793ba7696548616/cpp/dolfinx/geometry/utils.h#L543-L544

This call always uses `boost::multiprecision::cpp_bin_float_double_extended>` and in my experience is the time bottleneck inside of `determine_point_ownership`. Using `double` makes things much faster.

How would you go about exposing this template argument to users i.e. letting users decide what scalar type to use in `compute_distance_gjk`? Maybe an environment variable?

Here is a profiling example. This small script calls determine_point_ownership a few times:

```python3
from dolfinx import mesh, geometry
import numpy as np
from mpi4py import MPI
from line_profiler import LineProfiler

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

def main():
tdim = 3
num_elems_side = 20
domain = mesh.create_unit_cube(MPI.COMM_WORLD,
num_elems_side, num_elems_side, num_elems_side,
mesh.CellType.tetrahedron, dtype=np.float64)
call_determine_point_ownership(domain)

def call_determine_point_ownership(domain):
num_points = 100000
for it in range(10):
rng = np.random.default_rng(it)
points = rng.random((num_points, 3))
po = geometry._cpp.geometry.determine_point_ownership(
domain._cpp_object,
points,
1e-6)

if __name__=="__main__":
lp = LineProfiler()
lp.add_function(call_determine_point_ownership)
lp_wrapper = lp(main)
lp_wrapper()
with open(f"profiling_rank{rank}.txt", 'w') as pf:
lp.print_stats(stream=pf, sort=True, stripzeros=True)
```
With `boost::multiprecision` it takes 32.9 seconds, and with `double` instead, it takes about 3.1 seconds on my container. Results might be more accurate with `boost::multiprecision`, but doubles are more than enough for my particular application.

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.