Attempting to write a basic Hermitian eigenproblem

variational formulation
1

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

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
weakformimg743×62 6.16 KB

2

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.
3

Thanks for the reply. Here is the latest iteration

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)

4

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

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

    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.
1 Like