From 901ca699a407bead62454c30e28befdcb3a588dd Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 15 Dec 2025 15:12:31 +0000 Subject: [PATCH 1/5] add test --- .../firedrake/regression/test_interp_dual.py | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index 000b1023b2..d0cfbe5561 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -2,6 +2,7 @@ import numpy as np from firedrake import * from firedrake.utils import complex_mode +from firedrake.matrix import MatrixBase import ufl @@ -352,3 +353,47 @@ def test_interp_dual_mixed(source_space, target_space): assert result is tensor for x, y, in zip(result.subfunctions, expected.subfunctions): assert np.allclose(x.dat.data_ro, y.dat.data_ro) + + +def test_assemble_action_adjoint(V1, V2): + u = TrialFunction(V1) + + a = interpolate(u, V2) # V1 x V2^* -> R, equiv. V1 -> V2 + assert a.arguments() == (TestFunction(V2.dual()), TrialFunction(V1)) + + a_adj = adjoint(a) # V2^* x V1 -> R, equiv. V2^* -> V1^* + assert a_adj.arguments() == (TestFunction(V1), TrialFunction(V2.dual())) + + f = assemble(inner(1, TestFunction(V2)) * dx) + + expr = action(a_adj, f) + assert isinstance(expr, Action) + res = assemble(expr) + assert isinstance(res, Cofunction) + assert res.function_space() == V1.dual() + + expr2 = action(f, a) # This simplifies into an Interpolate + assert isinstance(expr2, Interpolate) + res2 = assemble(expr2) + assert isinstance(res2, Cofunction) + assert res2.function_space() == V1.dual() + assert np.allclose(res.dat.data, res2.dat.data) + + A = assemble(a) + assert isinstance(A, MatrixBase) + + # This explicitly assembles the adjoint of A, then matmults with f + expr3 = action(adjoint(A), f) + assert isinstance(expr3, Action) + res3 = assemble(expr3) + assert isinstance(res3, Cofunction) + assert res3.function_space() == V1.dual() + assert np.allclose(res.dat.data, res3.dat.data) + + # This doesn't explicitly assemble the adjoint of A, but uses multHermitian + expr4 = action(f, A) + assert isinstance(expr4, Action) + res4 = assemble(expr4) + assert isinstance(res4, Cofunction) + assert res4.function_space() == V1.dual() + assert np.allclose(res2.dat.data, res4.dat.data) From 2a4f0c2efd8e81255942acf808f0d56cf0b3e11a Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 15 Dec 2025 15:12:40 +0000 Subject: [PATCH 2/5] add assembly case update comment --- firedrake/assemble.py | 8 ++++++++ tests/firedrake/regression/test_interp_dual.py | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index ca5836be95..2c942c4f46 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -503,6 +503,14 @@ def base_form_assembly_visitor(self, expr, tensor, bcs, *args): with lhs.dat.vec_ro as x, rhs.dat.vec_ro as y: res = x.dot(y) return res + elif isinstance(rhs, matrix.MatrixBase): + # Compute action(Cofunc, Mat) => Mat^* @ Cofunc + petsc_mat = rhs.petscmat + (_, col) = rhs.arguments() + res = tensor if tensor else firedrake.Function(col.function_space().dual()) + with lhs.dat.vec_ro as v_vec, res.dat.vec as res_vec: + petsc_mat.multHermitian(v_vec, res_vec) + return res else: raise TypeError("Incompatible RHS for Action.") else: diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index d0cfbe5561..3584cf92e2 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -382,7 +382,7 @@ def test_assemble_action_adjoint(V1, V2): A = assemble(a) assert isinstance(A, MatrixBase) - # This explicitly assembles the adjoint of A, then matmults with f + # This (currently) explicitly assembles the adjoint of A, then matmults with f expr3 = action(adjoint(A), f) assert isinstance(expr3, Action) res3 = assemble(expr3) From 598bfe0055e13d6f8f86c0dfdf01dd7770766f79 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Mon, 15 Dec 2025 18:08:15 +0000 Subject: [PATCH 3/5] add simplification --- firedrake/assemble.py | 4 ++++ tests/firedrake/regression/test_interp_dual.py | 10 +++++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 2c942c4f46..c7d6a0d7c7 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -838,6 +838,10 @@ def restructure_base_form(expr, visited=None): # Replace arguments return ufl.replace(right, replace_map) + if isinstance(left, ufl.Adjoint) and isinstance(right, ufl.Cofunction) and left.arguments()[-1].function_space() == right.function_space(): + # Action(Adjoint(A), w*) -> Action(w*, A) + return ufl.action(right, left.form()) + # -- Case (4) -- # if isinstance(expr, ufl.Adjoint) and isinstance(expr.form(), ufl.core.base_form_operator.BaseFormOperator): B = expr.form() diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index 3584cf92e2..ee83d2b9e7 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -382,18 +382,18 @@ def test_assemble_action_adjoint(V1, V2): A = assemble(a) assert isinstance(A, MatrixBase) - # This (currently) explicitly assembles the adjoint of A, then matmults with f - expr3 = action(adjoint(A), f) + # This doesn't explicitly assemble the adjoint of A, but uses multHermitian + expr3 = action(f, A) assert isinstance(expr3, Action) res3 = assemble(expr3) assert isinstance(res3, Cofunction) assert res3.function_space() == V1.dual() assert np.allclose(res.dat.data, res3.dat.data) - # This doesn't explicitly assemble the adjoint of A, but uses multHermitian - expr4 = action(f, A) + # This is simplified into action(f, A) to avoid explicit assembly of adjoint(A) + expr4 = action(adjoint(A), f) assert isinstance(expr4, Action) res4 = assemble(expr4) assert isinstance(res4, Cofunction) assert res4.function_space() == V1.dual() - assert np.allclose(res2.dat.data, res4.dat.data) + assert np.allclose(res.dat.data, res4.dat.data) From 51b08e2217ed4c0778c04e6cae5aae965020990c Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 16 Dec 2025 11:05:04 +0000 Subject: [PATCH 4/5] fixes --- firedrake/assemble.py | 3 +- .../firedrake/regression/test_interp_dual.py | 70 +++++++++---------- 2 files changed, 36 insertions(+), 37 deletions(-) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index c7d6a0d7c7..14c737fe91 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -32,6 +32,7 @@ from pyop2.exceptions import MapValueError, SparsityFormatError from pyop2.types.mat import _GlobalMatPayload, _DatMatPayload from pyop2.utils import cached_property +from ufl.algorithms.analysis import extract_arguments __all__ = "assemble", @@ -838,7 +839,7 @@ def restructure_base_form(expr, visited=None): # Replace arguments return ufl.replace(right, replace_map) - if isinstance(left, ufl.Adjoint) and isinstance(right, ufl.Cofunction) and left.arguments()[-1].function_space() == right.function_space(): + if isinstance(left, ufl.Adjoint) and len(extract_arguments(right)) == 1: # Action(Adjoint(A), w*) -> Action(w*, A) return ufl.action(right, left.form()) diff --git a/tests/firedrake/regression/test_interp_dual.py b/tests/firedrake/regression/test_interp_dual.py index ee83d2b9e7..b0ce971a95 100644 --- a/tests/firedrake/regression/test_interp_dual.py +++ b/tests/firedrake/regression/test_interp_dual.py @@ -361,39 +361,37 @@ def test_assemble_action_adjoint(V1, V2): a = interpolate(u, V2) # V1 x V2^* -> R, equiv. V1 -> V2 assert a.arguments() == (TestFunction(V2.dual()), TrialFunction(V1)) - a_adj = adjoint(a) # V2^* x V1 -> R, equiv. V2^* -> V1^* - assert a_adj.arguments() == (TestFunction(V1), TrialFunction(V2.dual())) - - f = assemble(inner(1, TestFunction(V2)) * dx) - - expr = action(a_adj, f) - assert isinstance(expr, Action) - res = assemble(expr) - assert isinstance(res, Cofunction) - assert res.function_space() == V1.dual() - - expr2 = action(f, a) # This simplifies into an Interpolate - assert isinstance(expr2, Interpolate) - res2 = assemble(expr2) - assert isinstance(res2, Cofunction) - assert res2.function_space() == V1.dual() - assert np.allclose(res.dat.data, res2.dat.data) - - A = assemble(a) - assert isinstance(A, MatrixBase) - - # This doesn't explicitly assemble the adjoint of A, but uses multHermitian - expr3 = action(f, A) - assert isinstance(expr3, Action) - res3 = assemble(expr3) - assert isinstance(res3, Cofunction) - assert res3.function_space() == V1.dual() - assert np.allclose(res.dat.data, res3.dat.data) - - # This is simplified into action(f, A) to avoid explicit assembly of adjoint(A) - expr4 = action(adjoint(A), f) - assert isinstance(expr4, Action) - res4 = assemble(expr4) - assert isinstance(res4, Cofunction) - assert res4.function_space() == V1.dual() - assert np.allclose(res.dat.data, res4.dat.data) + f_form = inner(1, TestFunction(V2)) * dx + + for f in (f_form, assemble(f_form)): + expr = action(adjoint(assemble(a)), f) + assert isinstance(expr, Action) + res = assemble(expr) + assert isinstance(res, Cofunction) + assert res.function_space() == V1.dual() + + expr2 = action(f, a) # This simplifies into an Interpolate + assert isinstance(expr2, Interpolate) + res2 = assemble(expr2) + assert isinstance(res2, Cofunction) + assert res2.function_space() == V1.dual() + assert np.allclose(res.dat.data, res2.dat.data) + + A = assemble(a) + assert isinstance(A, MatrixBase) + + # This doesn't explicitly assemble the adjoint of A, but uses multHermitian + expr3 = action(f, A) + assert isinstance(expr3, Action) + res3 = assemble(expr3) + assert isinstance(res3, Cofunction) + assert res3.function_space() == V1.dual() + assert np.allclose(res.dat.data, res3.dat.data) + + # This is simplified into action(f, A) to avoid explicit assembly of adjoint(A) + expr4 = action(adjoint(A), f) + assert isinstance(expr4, Action) + res4 = assemble(expr4) + assert isinstance(res4, Cofunction) + assert res4.function_space() == V1.dual() + assert np.allclose(res.dat.data, res4.dat.data) From 8aa0267f3cf9a2ae9259ae421a51afacb4c2de95 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 16 Dec 2025 11:44:51 +0000 Subject: [PATCH 5/5] change check, add comment --- firedrake/assemble.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 14c737fe91..e495995897 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -32,7 +32,6 @@ from pyop2.exceptions import MapValueError, SparsityFormatError from pyop2.types.mat import _GlobalMatPayload, _DatMatPayload from pyop2.utils import cached_property -from ufl.algorithms.analysis import extract_arguments __all__ = "assemble", @@ -839,8 +838,10 @@ def restructure_base_form(expr, visited=None): # Replace arguments return ufl.replace(right, replace_map) - if isinstance(left, ufl.Adjoint) and len(extract_arguments(right)) == 1: - # Action(Adjoint(A), w*) -> Action(w*, A) + # Action(Adjoint(A), w*) -> Action(w*, A) + if isinstance(left, ufl.Adjoint) and not isinstance(right, firedrake.Function) and is_rank_1(right): + # TODO: ufl.action(Coefficient, Form) currently fails. When it is fixed, we can remove the + # `not isinstance(right, firedrake.Function)` check. return ufl.action(right, left.form()) # -- Case (4) -- #