Fix assembly and solves with the Hessian of an interpolation by danshapero · Pull Request #5478 · firedrakeproject/firedrake

danshapero · 2026-09-22T21:31:47Z

Resolves #5477: solving derivative(J, u) == 0 failed for a functional J containing interpolate(u, P) with P on a VertexOnlyMesh. The Hessian of the data term is Iᵀ M I.

Depends on FEniCS/ufl#532, which computes that Hessian symbolically. The DROP BEFORE MERGE commit installs that branch in CI.

Derivative expansion

preprocess_base_form expands the derivatives of a Form only if it has base form operators, whatever its mat_type. Plain Forms and Slate tensors are left to TSFC. Before, the "matfree" residual assembler skipped the expansion, so the chain rule through the interpolation produced a FormSum that OneFormAssembler rejects.
A "matfree" 2-form is still expanded only when its action is taken. MatrixFreeAssembler accepts any 2-form BaseForm, such as a FormSum, and keeps the tree intact so that a single implicit matrix assembles its action.

Dual slot of base form operators

UFL now simplifies Action(N(u; v̂), F) to N(u; F) for any base form operator and 1-form F, so the assembler can no longer split an operator from its dual slot by building an Action. Restructure case (5) and the Interpolate(v, 2-form) split are removed. Instead, an operator that cannot consume its assembled dual slot is assembled as N(u; v̂) and contracted with it:

an Interpolate whose dual slot assembles to a Matrix;
an ExternalOperator with a Cofunction or Matrix in its dual slot and no assembly method registered for it (assembly_method()).

The numerics of the Action branch move to contract, which both paths use.

Vertex-only meshes

Only the output matrix uses allocation_integral_types. An inner matrix allocated with the output's facet integral types crashed on a vertex-only mesh.

Tests

test_solve_interp_hessian: a Newton solve with the Hessian of an interpolation onto a vertex-only mesh.
test_matfree_form_sum_uses_implicit_matrix and test_matfree_slate_sum_product_uses_implicit_matrix.
test_reduced_functional.py::test_interpolate runs a second-order Taylor test with the Hessian.

🤖 Generated with Claude Code

Solving derivative(J, u) == 0 for a functional J that contains interpolate(u, P), with P a space on a vertex-only mesh, failed in two places in the base form assembler: - The residual assembler of the solver uses mat_type="matfree", which skipped the expansion of derivatives. The chain rule through the interpolation turns the residual into a FormSum, which the 1-form assembler rejects. The mat_type only applies to 2-forms, so derivatives of 0-forms and 1-forms are now always expanded. - Every inner matrix of the DAG was allocated with the integral types of the output matrix. The mass matrix on the vertex-only mesh then asked for interior facets, which that mesh does not have. Only the output matrix now uses those integral types. The test also needs FEniCS/ufl#525, which fixes the symbolic Hessian. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

pbrubeck · 2026-09-23T01:01:09Z

        original_expr = expr
-        if mat_type != "matfree":
-            # Don't expand derivatives if `mat_type` is 'matfree'
+        if mat_type != "matfree" or len(expr.arguments()) < 2:


don't we set mat_type="mafree" if len(expr.arguments()) < 2 already?

pbrubeck · 2026-09-23T01:03:09Z

preprocess_base_form now always expands the derivatives of 0-forms and 1-forms. The "matfree" mat_type of the residual assembler skipped this, so the chain rule through the interpolation gave a FormSum that OneFormAssembler rejects.

This is not ideal. We should only be expanding them as a last resort

pbrubeck · 2026-09-23T01:05:02Z

Only the output matrix uses allocation_integral_types. An inner matrix allocated with the output's facet integral types crashes on a vertex-only mesh.

I think the logic should involve checking for the vertex-only mesh, as opposed to changing the generic case

pbrubeck · 2026-09-23T01:07:35Z

Other tests seem to be breaking, please attribute failures and come back

pbrubeck · 2026-09-23T08:46:01Z

Only the output matrix uses allocation_integral_types. An inner matrix allocated with the output's facet integral types crashes on a vertex-only mesh.

I think the logic should involve checking for the vertex-only mesh, as opposed to changing the generic case

Actually this might be fine. We should not be propagating allocation_integral_types to the inner matrix

Install UFL from danshapero/fix-bfo-hessian (FEniCS/ufl#525) in the shared and documentation CI environments, so that CI tests this PR against the matching UFL fix. Drop this commit once firedrakeproject#525 merges. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…late Reconstructing an Interpolate now accepts the argument_slots keyword, and replaces a plain UFL Coargument in the dual slot with a Firedrake Argument on the dual of the target space. Also expand derivatives of matfree 0- and 1-forms that contain base form operators, so that the DAG traversal can assemble them. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…sian

pbrubeck · 2026-09-28T11:43:27Z

    assert np.allclose(u.dat.data, u2.dat.data)


+@pytest.mark.parallel([1, 3])


Suggested change

@pytest.mark.parallel([1, 3])

@pytest.mark.parallel([1, 3])

@pytest.mark.skipcomplex

connorjward · 2026-09-28T12:00:22Z

-            bcs = self._bcs if isinstance(e, ufl.BaseForm) and len(e.arguments()) == 2 else ()
-            return self.base_form_assembly_visitor(e, t, bcs, *operands)
+        if needs_matfree_assembler(self._form, self._mat_type, self._diagonal):
+            return self._matrix_free_assembler.assemble(tensor=tensor)


Can this check not happen inside the visitor in some way?

We need to check this early enough because this replaces a 2-form with its action (a 1-form). I wouldn't move it inside the visitor

Why not here? It's the next call site after this.

Actually, you are right. It makes a lot of sense to have this inside the visitor, to make sure inner matrices do not get explicitly assembled in matfree mode.

Putting this inside the visitor is harder, it will skip preprocessing steps.

How? It is effectively the same code just moved about.

The post-order visitor cannot safely be the dispatch point for MatrixFreeAssembler: it evaluates child operators before its callback, while an unexpanded matrix-free Hessian must reach the implicit matrix as a whole so its action can be preprocessed.

The right seam is the current early branch in BaseFormAssembler.assemble(), before _assemble_base_form; leaf two-forms already select their own assembler through TwoFormAssembler.

Also consider compositions like assemble(FormSum(Matrix, Matrix), mat_type="matfree"). It is much more simpler to apply the Action at the root, otherwise we would need a specialized traversal that applies the distributive property down to each two-from leaf.

Expanding the derivatives of a Form with base form operators may turn it into another BaseForm, which determines how it is assembled. Plain Forms and Slate tensors keep their type, so TSFC expands their derivatives. A 'matfree' matrix is still expanded when its action is taken. `preprocess_base_form` used `mat_type` to decide, so a 1-form with base form operators, whose default `mat_type` is 'matfree', was not expanded. `ufl.replace` then expanded it inside the assembly visitor, turning it into a FormSum. Expand it in preprocessing instead, and drop the second traversal in the visitor. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

FEniCS/ufl#525 was split into firedrakeproject#530, firedrakeproject#531 and firedrakeproject#532. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

pbrubeck · 2026-10-05T11:08:18Z

I pushed two commits here, to go with FEniCS/ufl#532 (#525 was split into #530, #531 and #532).

0f5df75: Only expand derivatives of forms with base form operators. In UFL, derivative() now only builds derivative nodes, and expand_derivatives is the only place where a derivative can turn a Form into another BaseForm. On this side, preprocess_base_form expands derivatives only where that can change how the form is assembled:

A Form with base form operators is expanded, whatever its rank.
Plain Forms and Slate tensors are left to TSFC, which also handles spatial derivatives.
A matfree 2-form is still expanded only when its action is taken.

Before, the check was mat_type != "matfree". A 1-form with base form operators defaults to matfree, so it was not expanded, and ufl.replace then expanded it inside the assembly visitor. That is why the visitor needed a second traversal, which this commit removes. With the new UFL, removing the second traversal without the new check makes test_solve_interp_hessian fail with a FormSum reaching OneFormAssembler.

7077538: DROP BEFORE MERGE. CI now installs FEniCS/ufl@pbrubeck/bfo-dual-slot-derivative instead of the closed danshapero/ufl@fix-bfo-hessian.

Tested locally with that UFL branch:

test_interp_dual, test_assemble_baseform, external_operators, test_interpolate_cross_mesh, test_reduced_functional and test_slate_infrastructure all pass, with nothing failing, on 1 and 3 ranks.
On 1 rank I also reran 22 more files from regression, slate and adjoint with and without this commit. The same tests failed on both sides. These failures come from my local build and FIAT branch, not from this change.

🤖 Generated with Claude Code

UFL now simplifies Action(N(u; v*), F) to N(u; F) for any base form operator N and 1-form F, so the assembler can no longer split an operator from its dual slot by building an Action: restructure case (5) and the Interpolate split folded straight back. Instead, when an operator cannot consume its assembled dual slot, assemble N(u; vhat) with a coargument vhat and contract it with the assembled dual slot: - An Interpolate whose dual slot assembles to a Matrix. - An ExternalOperator with a Cofunction or Matrix in its dual slot and no assembly method registered for it. `assembly_method()` looks this up. The numerics of the Action branch move to `contract`, which both paths use. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

An assembled Function has no arguments, so the contraction of an external operator with its assembled dual slot failed. Use the arguments of N(u; vhat) instead. Rename dual_slot_coargument to replace_dual_slot_by_coargument, which is what it does. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

danshapero added the LLM used An LLM was used in the production of this PR label Sep 22, 2026

pbrubeck reviewed Sep 23, 2026

View reviewed changes

danshapero and others added 11 commits September 23, 2026 06:54

Strengthen interpolation Taylor test

3446307

Use UFL backend-aware argument reconstruction

ac09410

BaseFormAssembler: fixes for matfree

446c184

Merge remote-tracking branch 'origin/main' into danshapero/interp-hes…

9b1080f

…sian

BaseFormAssembler: defer expansion for compilable trees

0b638b5

BaseFormAssembler: unify derivative expansion

14985a5

Continue derivative expansion after child preprocessing

9226394

revert expand_derivatives_form

b5efb0a

fixes

ca93e72

pbrubeck reviewed Sep 28, 2026

View reviewed changes

connorjward reviewed Sep 28, 2026

View reviewed changes

pbrubeck added 2 commits September 29, 2026 22:58

CI: retrigger Firedrake checks

8c7965f

Merge branch 'release' into danshapero/interp-hessian

8817aa2

pbrubeck changed the base branch from main to release October 5, 2026 09:21

pbrubeck added bug base:release Run this PR using a release build labels Oct 5, 2026

pbrubeck and others added 3 commits October 5, 2026 11:10

more tests

19b138e

DROP BEFORE MERGE: install UFL from FEniCS/ufl#532

7077538

FEniCS/ufl#525 was split into firedrakeproject#530, firedrakeproject#531 and firedrakeproject#532. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

pbrubeck mentioned this pull request Oct 6, 2026

Differentiate the dual slot of base form operators FEniCS/ufl#532

Open

		assert np.allclose(u.dat.data, u2.dat.data)


		@pytest.mark.parallel([1, 3])

	@pytest.mark.parallel([1, 3])
	@pytest.mark.parallel([1, 3])
	@pytest.mark.skipcomplex

Conversation

danshapero commented Sep 22, 2026 • edited by pbrubeck Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Derivative expansion

Dual slot of base form operators

Vertex-only meshes

Tests

Uh oh!

pbrubeck Sep 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pbrubeck commented Sep 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

pbrubeck commented Sep 23, 2026

Uh oh!

pbrubeck commented Sep 23, 2026

Uh oh!

pbrubeck commented Sep 23, 2026

Uh oh!

pbrubeck Sep 28, 2026

Choose a reason for hiding this comment

Uh oh!

connorjward Sep 28, 2026

Choose a reason for hiding this comment

Uh oh!

pbrubeck Sep 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

connorjward Sep 28, 2026

Choose a reason for hiding this comment

Uh oh!

pbrubeck Sep 28, 2026

Choose a reason for hiding this comment

Uh oh!

pbrubeck Sep 28, 2026

Choose a reason for hiding this comment

Uh oh!

connorjward Sep 28, 2026

Choose a reason for hiding this comment

Uh oh!

pbrubeck Oct 5, 2026

Choose a reason for hiding this comment

Uh oh!

pbrubeck commented Oct 5, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

danshapero commented Sep 22, 2026 •

edited by pbrubeck

Loading

pbrubeck Sep 23, 2026 •

edited

Loading

pbrubeck commented Sep 23, 2026 •

edited

Loading

pbrubeck Sep 28, 2026 •

edited

Loading