# Attempting to write a basic Hermitian eigenproblem

**URL:** <https://fenicsproject.discourse.group/t/attempting-to-write-a-basic-hermitian-eigenproblem/13191>\
**Category:** variational formulation\
**Created:** [January 1, 2024, 1:42am UTC](https://fenicsproject.discourse.group/t/attempting-to-write-a-basic-hermitian-eigenproblem/13191 "2024-01-01T01:42:27Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Morberticus\_Eab](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/morberticus_eab/32/2360_2.png) [@Morberticus\_Eab](https://fenicsproject.discourse.group/u/Morberticus_Eab)\
**Post date:** [January 1, 2024, 1:42am UTC](https://fenicsproject.discourse.group/t/attempting-to-write-a-basic-hermitian-eigenproblem/13191/1 "2024-01-01T01:42:27Z")

</div>

I’m trying to write a script to solve a basic Hermitian eigenproblem, ∇^2u = Eu, on a 2D square, using the new dolfinx module. I think I can follow the documentation once the problem is solved (extracting eigenvalues, plotting etc), but I am struggling with defining the problem. Here is what I have so far

```auto
from mpi4py import MPI                                                                         
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem, assemble_matrix
from ufl import ds, dx, grad, inner
 
msh = mesh.create_rectangle( 
    comm=MPI.COMM_WORLD,
    points=((0.0, 0.0), (1.0, 1.0)),
    n=(32, 32),
    cell_type=mesh.CellType.triangle
) 
V = fem.functionspace(msh, ("Lagrange", 1))
 
 # Set Dirichlet boundary conditions
def bc_checker(x):
    return np.logical_or(
        np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)),
        np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))
    )
 
facets = mesh.locate_entities_boundary(
    msh,
    dim=(msh.topology.dim - 1),
    marker=bc_checker
)
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
 
# Define problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(grad(u), grad(v)) * dx
 
A = assemble_matrix(a, bcs=[bc])
A.assemble()

solver = SLEPc.EPS().create(MPI.COMM_WORLD)
solver.setOperators(A)
solver.solve()
![img1|551x500](upload://sxJyQT10x2ZCqDX3CXPNUey2jWt.jpeg)

```

It fails at assemble\_matrix. “AttributeError: ‘Form’ object has no attribute ‘\_cpp\_object’”. I remember doing something similar above to hande the variational representation of the problem with the old dolfin module. Perhaps this is no longer the right procedure.

[edit] - For more info: This is the general form I am trying to reproduce

 ![weakformimg](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/7fe24c2cc72df33994af5f23dd511afb40adfab0.png)

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 1, 2024, 6:02am UTC](https://fenicsproject.discourse.group/t/attempting-to-write-a-basic-hermitian-eigenproblem/13191/2 "2024-01-01T06:02:10Z")

</div>

There are surely a couple of mistakes:

1. `A = assemble_matrix(a, bcs=[bc])` should be `A = assemble_matrix(fem.form(a), bcs=[bc])`
2. you need to assemble a matrix for the right-hand side and then solve a generalized eigenvalue problem, because the mass matrix is not the identity matrix for P1 finite elements.

---

<div class="post-metadata">

**Author:** ![Morberticus\_Eab](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/morberticus_eab/32/2360_2.png) [@Morberticus\_Eab](https://fenicsproject.discourse.group/u/Morberticus_Eab)\
**Post date:** [January 1, 2024, 2:27pm UTC](https://fenicsproject.discourse.group/t/attempting-to-write-a-basic-hermitian-eigenproblem/13191/3 "2024-01-01T14:27:59Z")

</div>

Thanks for the reply. Here is the latest iteration

```auto
from mpi4py import MPI
import numpy as np
import ufl
from slepc4py import SLEPc
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem, assemble_matrix
from ufl import ds, dx, grad, inner

msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=((0.0, 0.0), (1.0, 1.0)),
    n=(32, 32),
    cell_type=mesh.CellType.triangle
)
V = fem.functionspace(msh, ("Lagrange", 1))

 # Set Dirichlet boundary conditions
def bc_checker(x):
    return np.logical_or(
        np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)),
        np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))
    )

facets = mesh.locate_entities_boundary(
    msh,
    dim=(msh.topology.dim - 1),
    marker=bc_checker
)
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=0.0, dofs=dofs, V=V)

# Define problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(grad(u), grad(v)) * dx
m = u*v*dx

A = assemble_matrix(fem.form(a), bcs=[bc])
A.assemble()
M = assemble_matrix(fem.form(m), bcs=[bc])
M.assemble()

solver = SLEPc.EPS().create(MPI.COMM_WORLD)
solver.setOperators(A, M)
solver.solve()
nconv = solver.getConverged()

for i in range(nconv-1):
    print(solver.getEigenpair(i))

```

It solves the problem, but only for one eigenpair? I would like to solve for all eigenpairs (time efficiency is not an issue)

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [January 1, 2024, 4:44pm UTC](https://fenicsproject.discourse.group/t/attempting-to-write-a-basic-hermitian-eigenproblem/13191/4 "2024-01-01T16:44:54Z")

</div>

I have prepared a tutorial on this, since it was already on my todo list and it came up twice in two days. See [tutorial\_eigenvalue\_laplacian](https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_eigenvalue_laplacian.html)

You may want to limit yourself to `Approach #1` (and its fixup), since `Approach #2` requires a library you do not have (but could be easily installed, if needed).

If I extract the part in which `EPS` is setup

```auto
    eps = slepc4py.SLEPc.EPS().create(mesh.comm)
    eps.setOperators(A, B)
    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GHEP)
    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)
    eps.getST().getKSP().setType("preonly")
    eps.getST().getKSP().getPC().setType("lu")
    eps.getST().getKSP().getPC().setFactorSolverType("mumps")
    eps.solve()

```

you see:

- `eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GHEP)`: inform `SLEPc` that the problem is hermitian
- `eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)`: say that you want to compute 1 eigenvalue. If you need more, increase 1 to a larger number. However, be aware that when dealing with sparse matrices (as is always the case in FEM) you will typically not be able to compute all eigenvalues. I think that the default is computing 1 eigenvalue.
- `eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)`: say that the one eigenvalue you want is the smallest one, rather than the largest one. My understanding is that for this problem the eigenvalues of the problem in strong form go to +\infty, so there would be very little point in computing the largest one. I think that the default is computing the largest eigenvalue, so be careful to change this setting.