Constant initialization / assignment with division discards complex part

[jrmaddison]

e.g.

from firedrake import Constant
a = Constant(1.0 + 2.0j)
b = Constant(3.0 + 4.0j)
c = Constant(a / b)
d = Constant(0.0)
d.assign(a / b)
print(f"{complex(a) / complex(b)=}")
print(f"{complex(c)=} {complex(d)=}")

leads to

[...]/firedrake_complex/firedrake/src/ufl/ufl/algebra.py:257: ComplexWarning: Casting complex values to real discards the imaginary part
  e = float(a) / float(b)
complex(a) / complex(b)=(0.44+0.08j)
complex(c)=(0.3333333333333333+0j) complex(d)=(0.3333333333333333+0j)

Is this a UFL bug? Or should Firedrake lead to a TypeError being raised?

[wence-]

This is a UFL bug, but it's kind of subtle. UFL does the following to attempt to evaluate a numeric value for the division:

        a, b = self.ufl_operands
        a = a.evaluate(x, mapping, component, index_values)
        b = b.evaluate(x, mapping, component, index_values)
        # Avoiding integer division by casting to float
        try:
            e = float(a) / float(b)
        except TypeError:
            e = complex(a) / complex(b)
        return e

In this case a.evaluate(...) and b.evaluate(...) return scalar values of type numpy.complex128. While indeed float(1 + 2j) raises TypeError, float(numpy.complex128(1 + 2j)) returns 1 and emits a ComplexWarning.

Effectively none of the UFL scalar evaluation is type safe. The construction of c builds a symbolic a / b and then goes down this bad path. The assignment of a / b to d also does.

One way to fix this would be to override all the __dunder__ methods on firedrake Constant objects to handle the scalar division etc... cases by hand and then fall back to constructing symbolic expressions. But that's kind of just more sticking plaster.

[dham]

Could just fix UFL. Proactively check if either operand is complex instead of falling back.

[dham]

Actually, that integer division workaround is old Python 2 stuff. Probably the code should just call a/b