Eigenvalue solver

Question

Hi,

I'm using the SLEPc eigenvalue solver with a mass term of 300 (see below). When I use a coarse mesh, it produces a qualitatively incorrect eigenvector. When I change n (mesh points) to 100, it produces a correct eigenvector.

My question: Is there any option I can use to obtain a qualitatively correct eigenvector with a 10x10 mesh?

Thanks

from dolfin import *

n = 10

mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, 'CG', 1)
psi = Function(V)

def u0_boundary(x,on_boundary):
    return on_boundary
bc = DirichletBC(V,Constant(0.0),u0_boundary)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(Constant(300.0)*grad(u), grad(v))*dx + inner(u, v)*dx
m = inner(u, v)*dx
L = Constant(0.)*v*dx

A, M, _ = PETScMatrix(), PETScMatrix(), PETScVector()

assemble_system(a,L,bc,A_tensor=A,b_tensor=_)
assemble_system(m,L,A_tensor=M,b_tensor=_)

eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest real'
eigensolver.parameters['tolerance'] = 1.e-14
eigensolver.parameters['problem_type'] = 'gen_hermitian'

eigensolver.solve(5)

r, c, rx, cx = eigensolver.get_eigenpair(0)
print "%d:  %g" % (0, r)
psi.vector()[:] = rx

plot(psi)
interactive()

Answer 1

Hi, I set up the problem to look for eigenvalues closest to the smallest eigenvalue from the fine-grid solution.

from dolfin import *

n = 10

mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, 'CG', 1)
psi = Function(V)

def u0_boundary(x,on_boundary):
    return on_boundary
bc = DirichletBC(V,Constant(0.0),u0_boundary)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(Constant(300.)*grad(u), grad(v))*dx + inner(u, v)*dx
m = inner(u, v)*dx
L = Constant(0.)*v*dx

A, M, _ = PETScMatrix(), PETScMatrix(), PETScVector()

assemble_system(a,L,bc,A_tensor=A,b_tensor=_)
assemble_system(m,L,A_tensor=M,b_tensor=_)

eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'target real'
eigensolver.parameters['tolerance'] = 1.e-14
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 5924.   # Fine grid lambda_min

eigensolver.solve(5)

assert eigensolver.get_number_converged() > 0
r, c, rx, cx = eigensolver.get_eigenpair(0)
print "%d:  %g" % (0, r)
psi.vector()[:] = rx

plot(psi)
interactive()

The coarse-grid and fine-grid eigenvectors are quite similar. The error is mainly near the boundary.