firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
496 stars 157 forks source link

BUG: problem with .assign() #3448

Open hiroe- opened 5 months ago

hiroe- commented 5 months ago

I get ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() when doing

N = Function(W)
A = Function(W)
A.assign(weights[0]*N)

where W is a mixed function space:

V1 = FunctionSpace(mesh, "BDM", 2)
V2 = FunctionSpace(mesh, "DG", 1)
W = MixedFunctionSpace((V1, V2))

and weights is a list of float numbers. The full error message is below.

Traceback (most recent call last):
  File "/home/hyamazak/repos/sketchpad/averaging/disc_aver_sw_paradiag.py", line 1062, in <module>
    average(U0, V1, t=t)
  File "/home/hyamazak/repos/sketchpad/averaging/disc_aver_sw_paradiag.py", line 754, in average
    Average.assign(Average+dt*weights[0]*N)
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/hyamazak/firedrake/src/firedrake/firedrake/adjoint_utils/function.py", line 118, in wrapper
    ret = assign(self, other, *args, **kwargs)
  File "/home/hyamazak/firedrake/src/firedrake/firedrake/function.py", line 459, in assign
    elif expr == 0:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

The error disappears when using Constant() around weights[0]:

N = Function(W)
A = Function(W)
A.assign(Constant(weights[0])*N)

The relevant code is here at line 754:

https://github.com/colinjcotter/sketchpad/blob/a85a61b6b06c105d4f5aa3d45527893f94d910ce/averaging/disc_aver_sw_paradiag.py

To reproduce the error, get rid of Constant() from line 754 and run the code with

$ mpiexec -n 2 python disc_aver_sw_paradiag.py

Environment:

connorjward commented 5 months ago

So this is an interesting bit of Python and numpy magic that is going on. If you multiply a numpy float by a container (think list, tuple, array etc) then numpy actually goes and does the vectorised thing and multiplies by each element to produce a new array. E.g. np.float64(0.1) * [1, 2, 3] == [0.1, 0.2, 0.3].

I was able to produce a minimal example showing that custom data types need to define __len__ and __getitem__ for this to happen:

import numpy as np

class MyContainer:
    def __init__(self, data):
        self._data = data

    def __len__(self):
        return len(self._data)

    def __getitem__(self, index):
        return self._data[index]

    # Uncommenting this line prevents the behaviour
    # __iter__ = None

myobj = MyContainer([1, 2, 3])

print(np.float64(0.1) * myobj)  # prints [0.1, 0.2, 0.3]
print(0.1 * myobj)  # errors

What is happening specifically here is that Functions are iterable (i.e. you can write for indexed in my_function: ...) and so multiplying by a numpy data type unpacks this and an array is passed through to assign, which confuses it.

I think there are 3 solutions here:

  1. Document this as a confusing foot gun to avoid in assign.
  2. Make assign work with arrays.
  3. Make Functions not be iterable (set __iter__ to None).

I prefer option 3 but that would require changes in UFL that might be breaking. I don't think that it makes sense for Functions to be iterable.

@dham what is your opinion?