From b8c5252ab3bb02335ea3939c59f0e35e160fc610 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Wed, 18 Dec 2024 13:13:27 +0000 Subject: [PATCH 1/5] Fix for Cofunction self-assignment via interpolation --- firedrake/interpolation.py | 7 ++++++- tests/firedrake/regression/test_interp_dual.py | 7 +++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 0a598ba34b..449c4ff4d4 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -859,7 +859,12 @@ def _interpolate(self, *function, output=None, transpose=False, **kwargs): V = self.V result = output or firedrake.Function(V) with function.dat.vec_ro as x, result.dat.vec_wo as out: - mul(x, out) + if x.id != out.id: + mul(x, out) + else: + out_ = out.duplicate() + mul(x, out_) + out_.copy(result=out) return result else: diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index 5928d81e4e..ac24e3239f 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -34,6 +34,13 @@ def f1(mesh, V1): return Function(V1).interpolate(expr) +def test_interp_self(V1): + a = assemble(TestFunction(V1) * dx) + b = assemble(TestFunction(V1) * dx) + a.interpolate(a) + assert (a.dat.data_ro == b.dat.data_ro).all() + + def test_assemble_interp_operator(V2, f1): # Check type If1 = Interpolate(f1, V2) From 98f8c1ef5208a450b755133c8e6cfd4df7ac678a Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Wed, 18 Dec 2024 17:42:54 +0000 Subject: [PATCH 2/5] Update firedrake/interpolation.py Co-authored-by: David A. Ham --- firedrake/interpolation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 449c4ff4d4..04aced69d6 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -859,7 +859,7 @@ def _interpolate(self, *function, output=None, transpose=False, **kwargs): V = self.V result = output or firedrake.Function(V) with function.dat.vec_ro as x, result.dat.vec_wo as out: - if x.id != out.id: + if x is not out: mul(x, out) else: out_ = out.duplicate() From 63b74ca1bdcbb490e2d7f9b5db41da1c7a557a2c Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 19 Dec 2024 18:52:22 +0000 Subject: [PATCH 3/5] Add a more complicated test for dual interpolation where the tensor is a dependency of the input expression --- tests/firedrake/regression/test_interp_dual.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index ac24e3239f..e242951d74 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -41,6 +41,18 @@ def test_interp_self(V1): assert (a.dat.data_ro == b.dat.data_ro).all() +def test_assemble_interp_adjoint_tensor(mesh, V1, f1): + a = assemble(TestFunction(V1) * dx) + assemble(action(adjoint(Interpolate(f1 * TestFunction(V1), V1)), a), + tensor=a) + + x, y = SpatialCoordinate(mesh) + f2 = Function(V1, name="f2").interpolate( + exp(x) * y) + + assert np.allclose(assemble(a(f2)), assemble(Function(V1).interpolate(f1 * f2) * dx)) + + def test_assemble_interp_operator(V2, f1): # Check type If1 = Interpolate(f1, V2) From a2b6002a54d34b53f1e8986c44256c97ff29b546 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 19 Dec 2024 18:53:37 +0000 Subject: [PATCH 4/5] == -> np.allclose --- tests/firedrake/regression/test_interp_dual.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index e242951d74..73428e0cc6 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -38,7 +38,7 @@ def test_interp_self(V1): a = assemble(TestFunction(V1) * dx) b = assemble(TestFunction(V1) * dx) a.interpolate(a) - assert (a.dat.data_ro == b.dat.data_ro).all() + assert np.allclose(a.dat.data_ro, b.dat.data_ro) def test_assemble_interp_adjoint_tensor(mesh, V1, f1): From 2dc087aafb8d6d6ab0eba1b67f340a8e644a48eb Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 19 Dec 2024 18:57:23 +0000 Subject: [PATCH 5/5] Add an explanatory comment --- tests/firedrake/regression/test_interp_dual.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index 73428e0cc6..21c73c804c 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -43,6 +43,7 @@ def test_interp_self(V1): def test_assemble_interp_adjoint_tensor(mesh, V1, f1): a = assemble(TestFunction(V1) * dx) + # We want tensor to be a dependency of the input expression for this test assemble(action(adjoint(Interpolate(f1 * TestFunction(V1), V1)), a), tensor=a)