Eigenvalue solver problem

dolfinx
1

Hi everyone,
i have a question about defining the number of eigenvalues in a simple eigenvalue problem in dolfin-x.
in this code:

import numpy as np
import matplotlib.pyplot as plt
from dolfinx.plotting import plot
from mpi4py import MPI
from dolfinx import DirichletBC, Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
from dolfinx.cpp.mesh import CellType
import dolfinx
from dolfinx.io import XDMFFile
from math import sqrt
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (apply_lifting, assemble_matrix, assemble_vector,
                         locate_dofs_topological, set_bc)
nev=38

mesh = RectangleMesh(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([1, 1, 0])], [10, 10],
    cell_type=CellType.triangle)

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

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.setDimensions(nev)
eigensolver.setWhichEigenpairs(2)
eigensolver.setOperators(A)

eigensolver.solve()

for i in range (nev):
    l = eigensolver.getEigenpair(i)
    print(l.real)

there is some problems with some specific numbers.
for example for " nev = 38 " or " nev = 39 " there is this error:

---------------------------------------------------------------------------
Error                                     Traceback (most recent call last)
<ipython-input-10-8401c9868766> in <module>
     38 
     39 for i in range (nev):
---> 40     l = eigensolver.getEigenpair(i)
     41     print(l.real)

SLEPc/EPS.pyx in slepc4py.SLEPc.EPS.getEigenpair()

Error: error code 63
[0] EPSGetEigenpair() line 411 in /usr/local/slepc/src/eps/interface/epssolve.c
[0] Argument out of range
[0] Argument 2 out of range

but for the numbers lower than 38 or more than 40 it can compute the eigenvalues and there is not any error.
so, I can not understand the reason of the error for some specific numbers.
Thanks

2

Looks like you forgot to set the problem type. Adding

eigensolver.setProblemType(SLEPc.EPS.ProblemType.HEP)

will give some numbers at least.

1 Like
3

Dear @nate

Thanks for your help and sorry for the delay.

by adding Dirichlet Boundary Condition and two matrices as input parameters of “EPSSetOperators()”, the code can not runs properly when I am trying to collect the eigenvalues from smallest to largest values ( eigensolver.setWhichEigenpairs(2) ). But when I am trying to collect them from the largest to smallest values ( eigensolver.setWhichEigenpairs(1) ) there is not any problem and the code runs properly. I could not find any way to solve the problem.

In the following I have attached the codes and its outputs.

import numpy as np
import matplotlib.pyplot as plt
from dolfinx.plotting import plot
from mpi4py import MPI
from dolfinx import DirichletBC, Function, FunctionSpace, UnitSquareMesh, RectangleMesh
from dolfinx.fem import assemble_matrix, Form
from ufl import TrialFunction, TestFunction, grad, dx, dot
from slepc4py import SLEPc
from dolfinx.cpp.mesh import CellType
import dolfinx
from dolfinx.io import XDMFFile
from math import sqrt
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem import (apply_lifting, assemble_matrix, assemble_vector, locate_dofs_topological, set_bc)

nev=15

mesh = RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([0.01, 0.01, 0])], [100, 100], cell_type=CellType.triangle)

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

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

m0 = Function(V)

def boundary(x):
    return np.full(x.shape[1], True)
    
facets = locate_entities_boundary(mesh, 1, boundary)
bc = DirichletBC(m0, locate_dofs_topological(V, 1, facets))

A = assemble_matrix(a, bcs=[bc])
A.assemble()
B = assemble_matrix(b, bcs=[bc])
B.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev)
eigensolver.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eigensolver.setWhichEigenpairs(1)
eigensolver.setOperators(A,B)

eigensolver.solve()

for i in range (nev):
    l = eigensolver.getEigenvalue(i)
    print(l.real)

the result is:

12884606589.887491 
12884606589.860409 
12880572884.786907 
12880572884.766441 
12880160603.670452 
12880160602.680973 
12876235995.878403 
12876235994.537903 
12873321949.421896 
12873321949.389496 
12873206723.191608 
12873206720.327929 
12869827520.723207 
12869827520.722536 
12868483436.616455

​and by changing “eigensolver.setWhichEigenpairs(1)” to “eigensolver.setWhichEigenpairs(2)” to collect the eigenvalues from the smallest to largest values there is an error:

Error                                     Traceback (most recent call last)
<ipython-input-6-ae37569cdb16> in <module>
     47 
     48 for i in range (nev):
---> 49     l = eigensolver.getEigenvalue(i)
     50     print(l.real)

SLEPc/EPS.pyx in slepc4py.SLEPc.EPS.getEigenvalue()

Error: error code 63
[0] EPSGetEigenvalue() line 450 in /usr/local/slepc/src/eps/interface/epssolve.c
[0] Argument out of range
[0] Argument 2 out of range
4

I’m running into a similar problem right now. Are you sure that assembling both A & B with the same boundary conditions is right ? I don’t know what dolfinx does under the rug, but I have a feeling it could be wrong.

As for your current problem (maybe the term is not so appropriate now) I would use eigensolver.setWhichEigenpairs(eigensolver.Which.SMALLEST_MAGNITUDE) as well as calling eigensolver.getConverged() before trying to print anything.