From 06e01bf874978bb22d7d945ce003094e168e9c00 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Wed, 18 Dec 2024 13:13:27 +0000 Subject: [PATCH 1/9] 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 1686886f2a..6acfd4b6c4 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -877,7 +877,12 @@ def _interpolate(self, *function, output=None, transpose=None, adjoint=False, ** 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 b43fc99d9f..30a10ea790 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -35,6 +35,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 cec2433cf5c20f25e00b89c0c01b3f1162a37f78 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Wed, 18 Dec 2024 17:42:54 +0000 Subject: [PATCH 2/9] 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 6acfd4b6c4..b4b18d3cc1 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -877,7 +877,7 @@ def _interpolate(self, *function, output=None, transpose=None, adjoint=False, ** 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 2a581d4cf584f9ae49c7b87c70a71984f27812e6 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 19 Dec 2024 18:52:22 +0000 Subject: [PATCH 3/9] 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 30a10ea790..ee1e82cf2e 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -42,6 +42,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 d9ada401c9e51408656fcf11c1a6d7a8e57ff861 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 19 Dec 2024 18:53:37 +0000 Subject: [PATCH 4/9] == -> 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 ee1e82cf2e..0b97337a7b 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -39,7 +39,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 b9b277c21eea8b0a79f878f423ea11ae83dd4a56 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 19 Dec 2024 18:57:23 +0000 Subject: [PATCH 5/9] 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 0b97337a7b..b96db8eaff 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -44,6 +44,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) From 51e202b33afe3815f2d1069894d14ae6b6d18643 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Wed, 8 Jan 2025 17:06:33 +0000 Subject: [PATCH 6/9] Add missing conj needed for complex --- tests/firedrake/regression/test_interp_dual.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index b96db8eaff..9f40eb186b 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -36,14 +36,14 @@ def f1(mesh, V1): def test_interp_self(V1): - a = assemble(TestFunction(V1) * dx) - b = assemble(TestFunction(V1) * dx) + a = assemble(conj(TestFunction(V1)) * dx) + b = assemble(conj(TestFunction(V1)) * dx) a.interpolate(a) assert np.allclose(a.dat.data_ro, b.dat.data_ro) def test_assemble_interp_adjoint_tensor(mesh, V1, f1): - a = assemble(TestFunction(V1) * dx) + a = assemble(conj(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) From 853c4bb363086cd55ca0f3d50c23e0261181f9ee Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 9 Jan 2025 09:23:10 +0000 Subject: [PATCH 7/9] Add another missing conj --- 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 9f40eb186b..f2b52f2953 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -52,7 +52,7 @@ def test_assemble_interp_adjoint_tensor(mesh, V1, f1): f2 = Function(V1, name="f2").interpolate( exp(x) * y) - assert np.allclose(assemble(a(f2)), assemble(Function(V1).interpolate(f1 * f2) * dx)) + assert np.allclose(assemble(a(f2)), assemble(Function(V1).interpolate(conj(f1 * f2)) * dx)) def test_assemble_interp_operator(V2, f1): From 68673675865db36a377409eb7d2857a15f642792 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 23 Jan 2025 14:56:15 +0000 Subject: [PATCH 8/9] Work around Firedrake issue #3988 --- 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 f2b52f2953..5875230391 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -45,7 +45,7 @@ def test_interp_self(V1): def test_assemble_interp_adjoint_tensor(mesh, V1, f1): a = assemble(conj(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), + assemble(Interpolator(f1 * TestFunction(V1), V1).interpolate(a, adjoint=True), tensor=a) x, y = SpatialCoordinate(mesh) From 9c783cb731b3a5afa0daa9bf0849634d492236e6 Mon Sep 17 00:00:00 2001 From: "James R. Maddison" Date: Thu, 23 Jan 2025 14:56:39 +0000 Subject: [PATCH 9/9] Add extra assemble in test_assemble_interp_adjoint_complex --- 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 5875230391..b2d8f64463 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -126,7 +126,7 @@ def test_assemble_interp_adjoint_complex(mesh, V1, V2, f1): f1 = Constant(3 - 5.j) * f1 a = assemble(conj(TestFunction(V1)) * dx) - b = Interpolator(f1 * TestFunction(V2), V1).interpolate(a, adjoint=True) + b = assemble(Interpolator(f1 * TestFunction(V2), V1).interpolate(a, adjoint=True)) x, y = SpatialCoordinate(mesh) f2 = Function(V2, name="f2").interpolate(