Assemble system for eigenvalue problem

Question

Hi,

I am using slepc4py within FEniCS to solve a quadratic eigenvalue problem :

# a + b X + c X^2 where a, b, c are (symmetric) bilinear forms for Test and TrialFunctions
# bcs are Dirichlet boundary conditions
A = PETScMatrix()
B = PETScMatrix()
C = PETScMatrix()

A = assemble(a, tensor = A)
B = assemble(b, tensor = B)
C = assemble(c, tensor = C)

for bc in bcs:
 bc.apply(A)
 bc.apply(B)
 bc.apply(C)

 A_mat = A.mat()
 B_mat = B.mat()
 C_mat = C.mat()

 E = SLEPc.PEP().create()
 E.setOperators([A_mat,B_mat,C_mat])

I wonder whether the matrix A, B, and C remain symmetric after bc.apply() ? And if there is another solution to assemble the system with the boundary conditions (I thought using assemble_system(), but I can't figure out how because I get a different matrix format : <Matrix wrapper of <PETScMatrix of size ...>> instead of <PETScMatrix of size ...>, which I didn't manage to convert into petsc4py.PETSc.Mat?

Thank you in advance,

MiroK · Answer 1 · Mar 6, 2015

Hi, the matrices do not in general remain symmetric with assemble+bc.apply combination. The way to go is assemble_system. Your problem is due to the way you call the assembly routines.

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(3, 3)

V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

bc = DirichletBC(V, Constant(0), DomainBoundary())

a = inner(grad(u), grad(v))*dx
m = inner(u, v)*dx
L = inner(Constant(1), v)*dx

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

# See that assemble + bc.apply does not imply symmetry
assemble(a, A)
bc.apply(A)
print 'Is symmetric?', np.linalg.norm(A.array() - A.array().T) < 1E-10

# But assemble system does
assemble_system(a, L, bc, A_tensor=A, b_tensor=b)
assemble_system(m, L, bc, A_tensor=M, b_tensor=b)

print 'Is symmetric?', np.linalg.norm(A.array() - A.array().T) < 1E-10
print 'Is symmetric?', np.linalg.norm(M.array() - M.array().T) < 1E-10

from slepc4py import SLEPc

E = SLEPc.PEP().create()
E.setOperators([A.mat(), M.mat(), A.mat()])

Assemble system for eigenvalue problem

1 Answer