RF Eigenproblem #4529

Shawniz · 2025-08-29T15:11:36Z

This PR introduces a new eigenproblem/eigensolver abstraction that enables one to compute spectra of the linearisation of a sequence of Firedrake operations. The RFEigenproblem/RFEigensolver tandem closely mirror the LinearEigenproblem/LinearEigensolver tandem.
The main changes are introduced in the RFEigenproblem abstraction. Inspired by the interfacing and ability to define matrices implicitly introduced by @JHopeCollins in the JHopeCollins/tao-rf-mats branch, a new context RestrictedReducedFunctionalCtx object is used in the creation of the RestrictedReducedFunctionalMat which represents the linearisation of a sequence of FIredrake operations on a restricted function space.

connorjward

Some early feedback. I think @JHopeCollins needs to review this.

firedrake/rf_eigensolver.py

connorjward · 2025-09-01T09:03:48Z

firedrake/rf_eigensolver.py

+
+
+class RFEigenproblem:
+    """Generalised linear eigenvalue problem.


Docstring does not really describe the reduced functional aspect of this class.

I modified the docstring. Now I believe it is clearer.

firedrake/restricted_functional_ctx.py

connorjward · 2025-09-01T09:05:32Z

firedrake/restricted_functional_ctx.py

+    PETSc = None
+
+
+def new_restricted_control_variable(reduced_functional, function_space, dual=False):


Please use type hints throughout.

I added type hints where I believed arguments are unclear.

Thanks, but please include them everywhere that you've written a docstring.

firedrake/restricted_functional_ctx.py

firedrake/rf_eigensolver.py

connorjward · 2025-09-04T08:56:51Z

firedrake/rf_eigensolver.py

+
+
+class RFEigensolver:
+    r"""Solve a LinearEigenproblem.


Shouldn't this be a RFEigenProblem?

I think now it is clearer. Apologies, forgot to change.

connorjward · 2025-09-04T08:59:10Z

firedrake/rf_eigensolver.py

+
+
+class RFEigenproblem:
+    """Generalised eigenvalue problem.


I don't think that this line really describes what is going on. This is just copied from eigensolver.py.

Updated for more clarity.

firedrake/rf_eigensolver.py

connorjward · 2025-09-04T09:00:15Z

firedrake/rf_eigensolver.py

+    apply_riez : bool
+        If True, when defining A's adjoint action, we can output into the dual space.
+    action : str
+        Determine whether forward action is TLM or adjoint.


Please enumerate the valid options here.

The action docstring is now duplicated so this one can be removed.

JHopeCollins

This is really great, thanks for your contribution!

I've left a lot of comments/suggestions but the vast majority are minor style changes, in general I'm very happy with this.

If you have the time, then a little demo would be fantastic to have. It could just be showing the Helmholtz problem, but it will be much easier for people to find as a demo than in the tests.

tests/firedrake/adjoint/test_rf_eigensolver.py

JHopeCollins · 2025-09-04T12:05:46Z

tests/firedrake/adjoint/test_rf_eigensolver.py

+    problem_form = LinearEigenproblem(a, bcs=bc, restrict=True)
+    solver_form = LinearEigensolver(problem_form, n_evals=num_eigenvalues)
+    solver_form.solve()
+    nconv = solver_form.es.getConverged()


Should we be checking how many converged here?

I have tried but faced a bit of obstacles. The aim of the test is to check that two eigenvalue abstractions outputs the same results.

tests/firedrake/adjoint/test_rf_eigensolver.py

JHopeCollins · 2025-09-04T12:09:51Z

tests/firedrake/adjoint/test_rf_eigensolver.py

+    evals_rf = [
+        1 / solver_rf.eigenvalue(i) for i in range(0, min(num_eigenvalues, nconv))
+    ]


Why are we grabbing eigenvalues from the LinearEigenSolver but inverse eigenvalues from the RFEigenSolver?
Is this the right thing to do? If so, please add a comment explaining why.

Suggested change

evals_rf = [

1 / solver_rf.eigenvalue(i) for i in range(0, min(num_eigenvalues, nconv))

]

evals_rf = [

1 / solver_rf.eigenvalue(i) for i in range(min(num_eigenvalues, nconv))

]

This is done because the eigenvalue problem in LinearEigenproblem arises from the spatial discretisation of the bilinear forms arising from the Helmholtz PDE with forcing function f=cu for eigenvalue c. The solution operator in this case is a mapping from the forcing function to the solution. A wave-of-hands argument is that Jf=u with the f substituted from above, and using the fact that J is linear (which it is as the PDE is steady and linear) rearranges to `Ju = 1/cu'.

JHopeCollins · 2025-09-04T13:16:17Z

firedrake/rf_eigensolver.py

+        A = self._problem.A
+        M = self._problem.M
+        if not isinstance(A, PETSc.Mat):
+            raise TypeError(f"A must be a PETSc matrix, got {A.__class__}")
+        if not isinstance(M, PETSc.Mat) and (M is not None):
+            raise TypeError(f"M must be a PETSc matrix or None (identity), got {M.__class__}")
+
+        self.es.setDimensions(nev=self.n_evals, ncv=self.ncv, mpd=self.mpd)
+        self.es.setOperators(A, M)


The operators should be set in __init__, before the call to set_from_options. This is so that:

The solver actually has access to the Mats when everything is being set up during set_from_options,

We don't incorrectly/unnecessarily reset the operators every time we call solve.

Agreed, this change could probably be mirrored in eigensolver.py in the LinearEigensolver class.

Yes, if they are also set in solve in the LinearEigenSolver then that is also a bug.

They also need to be before the call to set_from_option in the init to make sure that it's all set up correctly.

JHopeCollins · 2025-09-04T13:22:24Z

firedrake/rf_eigensolver.py

+    if isinstance(variables, (tuple, list)):
+        if isinstance(variables[0], Function):
+            return tuple(Function(function_space).interpolate(v) for v in variables)
+        else:
+            return tuple(
+                Cofunction(function_space.dual()).interpolate(v) for v in variables
+            )
+    elif isinstance(variables, Function):
+        return Function(function_space).interpolate(variables)
+    else:
+        return Cofunction(function_space.dual()).interpolate(variables)


You can use pyadjoint.Enlist to handle whether we have a list or not. I think this is the right way to do it.

Suggested change

if isinstance(variables, (tuple, list)):

if isinstance(variables[0], Function):

return tuple(Function(function_space).interpolate(v) for v in variables)

else:

return tuple(

Cofunction(function_space.dual()).interpolate(v) for v in variables

)

elif isinstance(variables, Function):

return Function(function_space).interpolate(variables)

else:

return Cofunction(function_space.dual()).interpolate(variables)

from ufl.duals import is_primal

variables = Enlist(variables)

V = function_space if is_primal(variables[0]) else function_space.dual()

interp_vars = [Function(V).interpolate(v) for v in variables)

return variables.delist(interp_vars)

firedrake/rf_eigensolver.py

JHopeCollins · 2025-09-04T13:29:54Z

tests/firedrake/adjoint/test_rf_eigensolver.py

+    # Check adjoint action
+    problem_rf_true = RFEigenproblem(
+        J, bcs=bc, apply_riesz=True, restrict=restrict, action="adjoint"
+    )
+    problem_rf_false = RFEigenproblem(
+        J, bcs=bc, apply_riesz=False, restrict=restrict, action="adjoint"
+    )


Seeing as we're specifically testing the implicit matrix here, I think it'd be simpler to just construct them explicitly rather than going via the RFEigenProblem.

I actually thought that the code is less messy in this case. Otherwise, would have to do casework depending on restricted an unrestricted stuff. Also, would have to add functions to determine restricted space. Much more clunky whilst this happens under the hood of RFEigenproblem.

Co-authored-by: Josh Hope-Collins <[email protected]>

…ero more effective.

…redrake into rf_eigenproblem_clean

Co-authored-by: Josh Hope-Collins <[email protected]>

Shawniz added 8 commits August 29, 2025 17:52

Added rf_eigensolver to init.

70d6188

Restricted functional context.

091e5ff

RF eigensolver abstraction.

8258803

RF eigensolver tests.

7567ee4

Init so necessary new pyadjoint objects added.

dabb809

Corrected docstring in RRFC.

8ab669a

Remove wrong comment.

5a9e76f

Linted. Except wants me to remove things from init.

ffcdea3

connorjward requested changes Sep 1, 2025

View reviewed changes

JHopeCollins self-requested a review September 1, 2025 13:06

Shawniz added 7 commits September 4, 2025 09:11

Numpydoc followed. RFEigenproblem docstring made clearer.

949b9d7

Tidy up comments.

433d065

Stop inheriting OptionsManager.

7c365cf

Type-hint where needed. Added RFAction to __init__.

61f1208

Now added RFAction to init.

21c1983

Combined restricted_functional_ctx and rf_eigensolver.

c2f80ab

Removed the old import.

9343a04

connorjward requested changes Sep 4, 2025

View reviewed changes

Shawniz added 5 commits September 4, 2025 12:06

"=" changed to is.

30a91bf

Type-hinted necessary things.

ae46cdf

No double type-hinting.

d610892

Clearer docstring in RFEigensolver.

9ed9f14

Updated doc for more clarity.

58a40ad

JHopeCollins requested changes Sep 4, 2025

View reviewed changes

Shawniz and others added 6 commits September 4, 2025 19:12

Removed action and added options.

4578bc2

Ensure argument and docstring align.

88062bb

Initialise zero primal/dual space objects without interpolation.

2b4b899

Co-authored-by: Josh Hope-Collins <[email protected]>

Set operators in init so don't have to be reset in solve. Also init z…

ec9a731

…ero more effective.

Merge branch 'rf_eigenproblem_clean' of https://github.com/Shawniz/fi…

777d104

…redrake into rf_eigenproblem_clean

Using context manager when working with tape.

f63c1de

Co-authored-by: Josh Hope-Collins <[email protected]>



		class RFEigenproblem:
		"""Generalised linear eigenvalue problem.

		PETSc = None


		def new_restricted_control_variable(reduced_functional, function_space, dual=False):

RF Eigenproblem #4529

Are you sure you want to change the base?

RF Eigenproblem #4529

Uh oh!

Conversation

Shawniz commented Aug 29, 2025

Uh oh!

connorjward left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JHopeCollins left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!