Fix for Cofunction self-assignment via interpolation
#3939
Conversation
| a = assemble(TestFunction(V1) * dx) | ||
| b = assemble(TestFunction(V1) * dx) | ||
| a.interpolate(a) | ||
| assert (a.dat.data_ro == b.dat.data_ro).all() |
There was a problem hiding this comment.
Not sure if this should be a np.allclose test
There was a problem hiding this comment.
I think you're right as written. No flops should be performed so that should be an equality.
There was a problem hiding this comment.
I don't think it's actually no flops, it's interpolation which happens to be the identify except for roundoff.
There was a problem hiding this comment.
Changed to np.allclose, I'm pretty sure now that's correct
| a = assemble(TestFunction(V1) * dx) | ||
| b = assemble(TestFunction(V1) * dx) | ||
| a.interpolate(a) | ||
| assert (a.dat.data_ro == b.dat.data_ro).all() |
There was a problem hiding this comment.
I think you're right as written. No flops should be performed so that should be an equality.
Co-authored-by: David A. Ham <[email protected]>
| if x.id != out.id: | ||
| mul(x, out) | ||
| else: | ||
| out_ = out.duplicate() |
There was a problem hiding this comment.
Shouldn't a no-op be the right thing?
There was a problem hiding this comment.
I don't think so, e.g.
u = fd.Function(space, name="u")
v = fd.Function(space, name="v")
w = fd.Cofunction(space.dual(), name="w")
...
fd.assemble(fd.action(
fd.adjoint(fd.derivative(interpolate(v * u, space), u)),
w), tensor=w)
There was a problem hiding this comment.
I see, it only is a no-op in the MFE from the issue: assemble(L(w), tensor=w) with L being the identity.
There was a problem hiding this comment.
Maybe it'd be good to add a more complicated test then
There was a problem hiding this comment.
Done. I think there are also be issues with complex here (transpose should be Hermitian for the adjoint actions?).
…s a dependency of the input expression
Fixes #3935