Efficient way to convert data from numpy to fenics Function

1

Hi,

assuming a function

mesh = IntervalMesh(comm, 10, 0.0, 0.1)
V = FunctionSpace(mesh, 'CG', 1)
u_0 = Constant(1.0)
u_n = interpolate(u_0, V)

What would be the fastest way to convert its values to a numpy array and back?

Converting the function’s value to a numpy array:

In [5]: time u_n.vector()[:]
CPU times: user 85 µs, sys: 60 µs, total: 145 µs
Wall time: 160 µs
Out[5]: array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

time u_n.vector().get_local()
CPU times: user 58 µs, sys: 40 µs, total: 98 µs
Wall time: 105 µs
Out[6]: array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

In [7]: %timeit u_n.vector()[:]
4.85 µs ± 35.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [8]: %timeit u_n.vector().get_local()
5.01 µs ± 156 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Alternatively, using an Expression:

In [20]: expression_form = Expression("1.0", element=V.ufl_element())

n [21]: time interpolate(expression_form, V).vector()[:]
CPU times: user 528 µs, sys: 0 ns, total: 528 µs
Wall time: 537 µs
Out[21]: array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

In [22]: time interpolate(expression_form, V).vector().get_local()
CPU times: user 521 µs, sys: 0 ns, total: 521 µs
Wall time: 530 µs
Out[22]: array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

In [23]: %timeit interpolate(expression_form, V).vector()[:]
97.9 µs ± 1.92 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [24]: %timeit interpolate(expression_form, V).vector().get_local()
95.1 µs ± 854 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Defining the value of the function from a numpy array:

In [10]: %timeit  u_n.vector()[:] = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
36.9 µs ± 1 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [12]: %timeit  u_n.vector().get_local()[:] = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
8.21 µs ± 13.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

As my model demands several of these convertions, I’m wondering if there’s a more efficient way of doing it to improve performance. Any idea/thought about this is appreciated. Thanks!

2

If you’re this concerned about performance, I’d recommend manipulating the data using PETSc directly, and avoiding python altogether.

3

This has been discussed in the fenics slack channel:

  1. set_local/vector the problem is set_local expects List[float] and when you pass F this F has magic getitem python method. It is implicitly lazy-calling F[0], F[1], F[2], … to get that List[float] and this is passed to set_local. Not only is List[float] not local in memory but its terrible to iterate over etc etc.
    (You can call set_local on some wrong arguments to see all possible signatures)
    set_local: 1. (self: dolfin.cpp.la.GenericVector, arg0: List[float]) → None

  2. set_local/get_local is a bit better, because F.get_local will be numpy array (yupii), however it will copy underlying data array to produce that numpy array. It is copy operation and also probably a conversion of numpy array to List[float] what you see here.

  3. array/get_local even better, because vector setitem interface allows to be set from numpy array, the only thing this you see here is copying numpy array on F.get_local call.
    setitem: 3. (self: dolfin.cpp.la.GenericVector, arg0: slice, arg1: numpy.ndarray[float64]) → None

  4. array/vector even better, because setitem has also following signature
    setitem: 2. (self: dolfin.cpp.la.GenericVector, arg0: slice, arg1: dolfin.cpp.la.GenericVector) → None
    which means you are assigning a Vector directly, i.e. doing just one copy assignment operation in cpp layer.

4

Thanks, but avoiding python is not an option for me …

5

Thanks a lot for your answer!