Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Constant initialization / assignment with division discards complex part #2376

Open
jrmaddison opened this issue Mar 7, 2022 · 4 comments · May be fixed by #4010
Open

Constant initialization / assignment with division discards complex part #2376

jrmaddison opened this issue Mar 7, 2022 · 4 comments · May be fixed by #4010

Comments

@jrmaddison
Copy link
Contributor

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-
Copy link
Contributor

wence- commented Mar 16, 2022

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
Copy link
Member

dham commented Mar 16, 2022

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

@dham
Copy link
Member

dham commented Mar 16, 2022

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

@jrmaddison
Copy link
Contributor Author

Fixed by FEniCS/ufl#342

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants