Symmetry and Dirichlet Eigenvalue problem

Question

I'm considering the Dirichlet Laplacian eigenvalue problem on some reasonable domain, and there are some things I'm confused about. Suppose I consider the unit square in 2D, for which I have closed form solutions. In that case, the following code gives me the answer that I expect:

mesh = UnitSquareMesh(10,10)

V = FunctionSpace(mesh, "Lagrange", degree = 1)
bc = DirichletBC(V, 0.0, 'on_boundary')
u = TrialFunction(V)
v = TestFunction(V)

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

L = inner(Constant(1), v)*dx

A = PETScMatrix()
assemble(a, A)
bc.apply(A)

M = PETScMatrix()
assemble(m, M)

#  Perform SLEPc eigenvalue solve with stiffness matrix A and mass matrix M\

However, the matrix A, after applying the Dirichlet boundary conditions, is no longer symmetric, which puzzles me. I would think that we should be making every effort to preserve such structure.

Answer 1

You should use assemble_system, which preserves symmetry.

Comments

When I use assemble_system, as was suggested elsewhere, I get zero pivot errors when I try to solve the eigenvalue problem.

---

Hi, perhaps you are not using the function correctly. The following works fine

from dolfin import *
from scipy.linalg import eigvalsh
from numpy import all as np_all

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, Constant(0), 'on_boundary')
u = TrialFunction(V)
v = TestFunction(V)

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

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

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

assert np_all(eigvalsh(A.array(), M.array())) 

---

That returns a bunch of eigenvalues with value = 1. The eigenvalues of this problem are of the form pi^ * (n^2 + m^2)

---

Remove the bcs on M or have a look at the eigenvalues larger than one.

---

Yes, this seems to be a resolution. If I use

assemble(m, M)

I get the correct eigenvalues and all of my matrices are symmetric.

Answer 2

from dolfin import *                                                            
import scipy.sparse as sp                                                       
from scipy.sparse.linalg import eigs                                            

mesh = UnitSquareMesh(10,10)                                                    

V = FunctionSpace(mesh, "Lagrange", degree = 1)                                 
bc = DirichletBC(V, 0.0, 'on_boundary')                                         
u = TrialFunction(V)                                                            
v = TestFunction(V)                                                             

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

L = inner(Constant(1), v)*dx                                                    

A = PETScMatrix()                                                               
assemble(a, tensor=A)                                                           
bc.apply(A)                                                                     

M = PETScMatrix()                                                               
assemble(m, tensor=M)                                                           

# This just converts PETSc to CSR                                               
A = sp.csr_matrix(A.mat().getValuesCSR()[::-1])                                 
M = sp.csr_matrix(M.mat().getValuesCSR()[::-1])                                 

v, V = eigs(A, 6, M, sigma=1.5)                                                 

print v

This calculates the physical eigenvalues of your problem. If you want to apply boundary conditions to your mass matrix M you can do a little trick to avoid the 1.0 eigenvalues in your spectrum.

from dolfin import *                                                                           
import scipy.sparse as sp                                                                      
import numpy as np                                                                             
from scipy.sparse.linalg import eigs                                                           

mesh = UnitSquareMesh(16,16)                                                                   

V = FunctionSpace(mesh, "Lagrange", degree = 1)                                                
bc = DirichletBC(V, 0.0, 'on_boundary')                                                        
u = TrialFunction(V)                                                                           
v = TestFunction(V)                                                                            

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

L = inner(Constant(1), v)*dx                                                                   

A = PETScMatrix()                                                                              
assemble(a, tensor=A)                                                                          
bc.apply(A)                                                                                    

M = PETScMatrix()                                                                              
assemble(m, tensor=M)                                                                          
# BCs applied this time                                                                        
bc.apply(M)                                                                                    

# If you had multiple bcs                                                                      
bcs = [bc]                                                                                     

bcinds = []                                                                                    
for bc in bcs:                                                                                 
    bcdict = bc.get_boundary_values()                                                          
    bcinds.extend(bcdict.keys())                                                               


# This just converts PETSc to CSR                                                              
A = sp.csr_matrix(A.mat().getValuesCSR()[::-1])                                                
M = sp.csr_matrix(M.mat().getValuesCSR()[::-1])                                                

# Create shift matrix                                                                          
shift = 1.2345e10*np.ones(len(bcinds))                                                         
S = sp.csr_matrix((shift, (bcinds, bcinds)), shape=A.shape)                                    

v, V = eigs(A+S, 6, M, sigma=1.0)                                                              

print v        

This will shift all your eigenvalues from arising from dirichlet boundary condition application to 1.2345e10. You can vary that value depending on which ones you need.