Add high-level interface to `SLEPc EPS`

[francesco-ballarin]

Describe new/missing feature

We now have a dolfinx.fem.petsc.LinearProblem and dolfinx.fem.petsc.NonlinearProblem.
It would be nice to also have a dolfinx.fem.slepc.LinearEigenvalueProblem that interfaces EPS.

I will work on this in the next few days.

Suggested user interface

[jhale]

How far down the rabbit hole do we want to go with adding interfaces? LinearProblem and NonlinearProblem were added because the first is common place for beginners and the second non-trivial. EPS seems, on the face of it, to be neither.

[francesco-ballarin]

How far down the rabbit hole do we want to go with adding interfaces?

That's a reasonable question. Surely, at least at this stage, not too far down the rabbit hole, for instance with TS or TAO. I'd think EPS is not as far down the rabbit hole as TS/TAO, but I'd like to hear more thoughts about this.

[jhale]

TAO is more niche. I'd imagine that writing one's own assemblers is non-trivial, so we should see how the recently added assemble_jacobian and assemble_residual work there and if anything can be modified. @schnellerhase can probably advise.

I've never used TS so I couldn't comment.

[jorgensd]

TAO is more niche. I'd imagine that writing one's own assemblers is non-trivial, so we should see how the recently added assemble_jacobian and assemble_residual work there and if anything can be modified. @schnellerhase can probably advise.

I've never used TS so I couldn't comment.

I tried to implement som TS stuff prior to the conference, and it was trickier than I hoped for.

I think having an EPS interface makes sense, as eigenvalue problems are commonplace.

[schnellerhase]

TAO assemblers are one to one the SNES ones, at least for the unconstrained/bound-constrained setting. Only difference is the additional objective that needs to be provided. For general constrained problems this changes and it requires quite some additional logic on top of 'just assembling'. Some of the specialised solvers from TAO also rely on different callbacks or callback interpretations, which makes a uniform interface to those tricky.

[michalhabera]

Please also have a look at https://github.com/fenics-dolfiny/dolfiny/blob/main/src/dolfiny/slepcproblem.py

Three points to consider:

1. It is generalized eigenvalue problems $T(\lambda)x = 0$, that are of real importance for engineers. SLEPc provides these, so one should not restrict from using these.
2. I find it more challenging to design how to interface nicely these from UFL. What we do in dolfiny we use trick with a placeholder function for $\lambda$, say lmbda0 = Function(...) so user provides a "residual form" for $T$ with lmbda0 possibly present there. Then for $T(\lambda) = A + \lambda B$ we collect the constant and linear terms as $A = T\rvert_{\lambda=0}$ and $B = \frac{\partial T}{\partial \lambda}\rvert_{\lambda=0}$. This works quite nicely as lmbda0 could be thought of as a parameter of the system for which the stability should be explored. I think @cmaurini has more experience here.
3. Dirichlet boundary conditions need special handling to not pollute the spectrum with BCs eigvals and modes. Firedrake added proper restriction approach recently I think, for this very reason.

[francesco-ballarin]

It is generalized eigenvalue problems that are of real importance for engineers. SLEPc provides these, so one should not restrict from using these.

I understand the point, but I wouldn't plan on doing it, I would only plan for a linear eigenvalue solver for now.

Dirichlet boundary conditions need special handling to not pollute the spectrum with BCs eigvals and modes. Firedrake added proper restriction approach recently I think, for this very reason.

That's the reason I was originally interested in this. In multiphenicsx we have a tutorial that uses restriction to throw away Dirichlet BCs dofs. I think this is one instance in which the multiphenicsx/restriction way of doing things and the submesh/mixed-dimensional approach are not functionally equivalent. DOFs restrictions wouldn't be incorporated, but at the very least there would be a Note: informing the user about spurious eigenpairs introduced by BCs.

[PretTy-Fish]

In multiphenicsx we have a tutorial that uses restriction to throw away Dirichlet BCs dofs.

I agree with this method and independently developed it while working on my own project before. It is more efficient as it effectively reduces the matrix size by removing the trivial block. It is unfortunate to hear that it wouldn't be incorporated, but may I ask where I can learn about the "submesh/mixed-dimensional approach"?

but at the very least there would be a Note: informing the user about spurious eigenpairs introduced by BCs.

I think it would also be useful to include the instruction to avoid them by changing the diagonal values of $A$ and $B$ when assembling, like in your multiphenicsx example and Jørgen's dolfinx_mpc tutorial. In the latter, it says:

The generalized eigenvalue problem will have spurious eigenvalues at lambda_spurious = diagval_A/diagval_B.

In personal experience, I found setting diagval_A = 1 and diagval_B = 0 moves the spurious eigenvalues to infinity, which works well when using the shift-invert method and does not mess up the eps.computeError() result too much. Setting them the other way around sets the spurious eigenvalue to 0, which helps when not using the shift-invert method.