A simple eigenvalue solver in DOLFIN-X

1

Hi,
I am trying to write a simple eigenvalue solver code in dolfin-x, but I have a problem in how to define a matrix as a PETScMatrix() and assembling the form into it. There is a proposed code in dolfin which I am trying to write it in dolfin-x:

from dolfin import *
mesh = …….
V = FunctionSpace(mesh, "Lagrange", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

# Assemble stiffness form
A = PETScMatrix()                                                      
assemble(a, tensor=A)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A)

# Compute all eigenvalues
eigensolver.solve()

# Extract largest (first) eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(0)

print "Largest eigenvalue: ", r

in the following, There is an ambiguity how to # create eigensolver and # compute the eigenvalues in DOLFIN-X? Thanks in advance.

1 Like
2

Hi Karo,

You can assemble the PETSc matrix by using the dolfinx.fem.assemble_matrix function in dolfinx. An example is given in this demo (https://github.com/FEniCS/dolfinx/blob/master/python/demo/stokes-taylor-hood/demo_stokes-taylor-hood.py). Refer to line 435-436.

Once you get the PETSc matrix. You can use slepc4py to compute the eigenvalues. You can refer to this slepc4py demo (https://gitlab.com/slepc/slepc4py/-/blob/master/demo/ex1.py) and adapt it to your need.

3 Likes
3

Disclaimer: This is something I tested quite sometime back with dolfinx-2019.2.9.99 so no guarantees that will work as is, but should be a good starting point in case if it doesn’t too

import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh
from dolfinx.fem import assemble_matrix, Form
from dolfinx.plotting import plot
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
# form = Form(a)

# Assemble stiffness tensor
A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)

eigensolver.solve()
vr, vi = A.getVecs()

lmbda = eigensolver.getEigenpair(0, vr, vi)
u = Function(V)
u.vector.setArray(vr.array)

# Plot eigenfunction (0-th)
plt.figure(figsize=(8,8))
eig1 = plot(u, cmap=plt.cm.jet)
plot(mesh)
plt.colorbar(eig1)

dolfinxNew
dolfinxNew800×800 78.9 KB

translated from

import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)

# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)
eigensolver = SLEPcEigenSolver(A)
eigensolver.solve()
r, c, rx, cx = eigensolver.get_eigenpair(0)
u = Function(V)
u.vector()[:] = rx

# Plot eigenfunction
plt.figure(figsize=(8,8))
eig1 = plot(u, cmap=plt.cm.jet)
plot(mesh)
plt.colorbar(eig1)

dolfinOld
dolfinOld800×800 78.7 KB

4 Likes
4

Dear @adeebkor
I read mentioned demos, thanks for your kind guidance.

5

Dear @bhaveshshrimali
Thanks for your help, it works for me properly too, and it would be really helpful.

6

Dear @bhaveshshrimali, @adeebkor
As I found, vr, vi = A.getVecs() extracts the real and complex part of the corresponding eigenvector. Do you have an idea how is it possible to extract the real and complex part of the eigenvalue in dolfin-x?
In dolfin, r, c, rx, cx = eigensolver.get_eigenpair(0), extracts both eigenvalue and eigenvector components.

7

lmbda = eigensolver.getEigenpair(0, vr, vi) and lmbda.real, lmbda.imag gives the real and imaginary part of eigen values. Check the example above for using slepc

8

Dear Karo,

vr, vi = A.getVecs() creates a placeholder vector to store the eigenvector of the matrix A. Following on the comment by bhaveshshrimali, once you called lmbda = eigensolver.getEigenpair(0, vr, vi), it will gives you the first eigenvalue and the first eigenvector. If you would like the rest eigenvalue, you can use loop just like the slepc4py example above.

Here is the documentation on slepc4py (https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS-class.html#getEigenpair) for the getEigenpair method.

9

Dear @adeebkor, @bhaveshshrimali,
Thanks again for your kind reply.
It is exactly my new problem about extracting more eigenvalues.
In this code, dolfin returns selected numbers of eigenvalues:

import numpy as np
import matplotlib.pyplot as plt
from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = PETScMatrix()
assemble(a, tensor=A)

eigensolver = SLEPcEigenSolver(A)
eigensolver.solve()

for i in range (10):
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print(r)

The result is:

7.822379233170962
7.5635582373483015
7.5635582373482855
7.306318516404323
7.15614193137886
7.156030031719445
6.901377998228639
6.901377998228612
6.6339159619522885
6.6339159619522565

But dolfin-x just returns the first eigenvalue:

import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
import dolfinx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)

eigensolver.solve()
self = eigensolver.getConverged()
vr, vi = A.createVecs()
print( "Number of converged eigenpairs %d" % self )
if self > 0:
    for i in range (self):
            l = eigensolver.getEigenpair(i ,vr, vi)
            print(l.real)

and The result is:

Number of converged eigenpairs 1
7.822379233170955

and it doesn’t return more eigenvalues. Please let me know if there is a solution to solve this problem. Thanks in advance.

10

You can set the number of eigenvalues that you want using the setDimensions() method. Here is your code that I’ve modified:

import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
import dolfinx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.setDimensions(nev=10)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = A.createVecs()
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    for i in range (evs):
            l = eigensolver.getEigenpair(i ,vr, vi)
            print(l.real)

You can refer to the documentation for more info.

2 Likes
11

Thanks a lot dear @adeebkor, @bhaveshshrimali,
The last topic about eigenvalue solver which I am interested in is :

esolver.parameters["spectrum"] = "smallest real"

this command helps to arrange eigenvalues from the smallest real part.
But I couldn’t find the command on the subject for dolfin-x, so please indicate if you have a source for this.

12

Hi,
A vague guess would be to just extract the eigen values using esolver.getEigenvalue(k) which should extract the k eigen values and then sort them yourself. However, looking up the corresponding function in old dolfin documentation for what it is doing under the hood should give you an idea. See this for the corresponding method, which should roughly translate to

eigensolver.setWhichEigenpairs(eigensolver.Which.LARGEST_REAL)
eigensolver.solve()

Adeeb can maybe provide a more detailed reply

1 Like
13

Most of the commands that I’ve wrote are slepc4py command. You need to look in the slepc4py documentation which is a separate package from dolfin-x.

Nevertheless, the command that you are looking for is eigensolver.setWhichEigenpairs(). This is a modified code that may return the smallest real eigenvalue that you need.

import numpy as np
from mpi4py import MPI
from dolfinx import Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
import dolfinx

mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = assemble_matrix(a, [])
A.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.setDimensions(nev=10)
eigensolver.setWhichEigenpairs(4)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = A.createVecs()
print( "Number of converged eigenpairs %d" % evs )
if evs > 0:
    for i in range (evs):
            l = eigensolver.getEigenpair(i ,vr, vi)
            print(l.real)

This page list the desired piece of spectrum (https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS.Which-class.html) and this page (https://slepc.upv.es/slepc4py-current/docs/apiref/slepc4py.SLEPc.EPS-class.html) is the documentation for the SLEPc.EPS class.

14

Dear @adeebkor, @bhaveshshrimali,
Thank you both for your kind help.