Eigenvalues of the Laplacian on the Sphere

Legacy DOLFIN
1

Hi all,

I am attempting to use DOLFIN to approximate the eigenvalues of the Laplacian on a sphere. I am trying to do this by calculating the eigenvalues on a 2pi x pi rectangle, with periodic boundary conditions on the left and right boundaries, and “constant” boundary conditions on the top and bottom boundaries (i.e., the entire bottom boundary needs to take on the same value. Same for the top). However, I seem to be getting imaginary eigenvalues when I run my code. Any ideas?


from dolfin import *
from dolfin import PETScMatrix, as_backend_type
from petsc4py import PETSc
import scipy.sparse as sp
from scipy.sparse.linalg import eigs
import numpy as np

# The rectangle will be divided into DIMENIONS x DIMENIONS regions
DIMENSIONS = 16

# Periodic BCs on the left + right boundaries
class PeriodicBoundary(SubDomain):

    def inside(self, x, on_boundary):
        return bool(near(x[0], 0) and on_boundary)
    
    def map(self, x, y):
        y[0] = x[0] - 2*DOLFIN_PI
        y[1] = x[1]
pbc = PeriodicBoundary()

# Create the mesh (2pi x pi rectangle) and the function space (Lagrange elements) constrained to the periodic BCs defined earlier
mesh = RectangleMesh(Point(0.0, 0.0, 0.0), Point(2*DOLFIN_PI, DOLFIN_PI, 0.0), DIMENSIONS, DIMENSIONS)
V = FunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

# Create our functions
u = TrialFunction(V)
v = TestFunction(V)

# Create our bilinear forms
a = dot(grad(u), grad(v))*dx
b = u*v*dx

# Get the matrix A and M corresponding to our bilinear forms
A = PETScMatrix()
M = PETScMatrix()
assemble(a, tensor=A)
assemble(b, tensor=M)

# Get the degrees of freedom (i.e., coordinates of vertices)
x = V.tabulate_dof_coordinates()

# Define MAT_LEN to be the amount of degrees of freedom
# Or, in other words, the number such that the matrices A and M have dimension MAT_LEN x MAT_LEN
MAT_LEN = 0
for i in x:
    MAT_LEN += 1

# bottom_elements - Used to store the degree of freedom numbers corresponding to those lying on the bottom side of the rectangle
# mat_index - The degree of freedom number (corresponds to the element of the array x)
# element_index - Used to index bottom_elements
bottom_elements = np.zeros(DIMENSIONS)
mat_index = 0
element_index = 0
# Loop through the degrees of freedom coordinates
for i in x:
    # If the y-coordinate of the degree of freedom is 0 (i.e. lying on the bottom side of the rectangle),
    #     then put that degree of freedom number (mat_index) in the array bottom_elements at the element_index-th spot. 
    #     Should loop DIMENSIONS times
    if (i[1] == 0):
        bottom_elements[element_index] = mat_index
        element_index += 1
    mat_index += 1

# top_elements - Used to store the degree of freedom numbers corresponding to those lying on the top side of the rectangle
# mat_index - The degree of freedom number (corresponds to the element of the array x)
# element_index - Used to index top_elements
top_elements = np.zeros(DIMENSIONS)
mat_index = 0
element_index = 0
# Loop through the degrees of freedom coordinates
for i in x:
    # If the y-coordinate of the degree of freedom is pi (i.e. lying on the top side of the rectangle),
    #     then put that degree of freedom number (mat_index) in the array top_elements at the element_index-th spot. 
    #     Should loop DIMENSIONS times
    if (i[1] == DOLFIN_PI):
        top_elements[element_index] = mat_index
        element_index += 1
    mat_index += 1

# Make a petsc matrix from the fenics matrices A and M. This allows us to make changes to the matrices using the petsc interface
A_petsc = as_backend_type(A).mat()
A_petsc.assemble()
A_petsc.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)

M_petsc = as_backend_type(M).mat()
M_petsc.assemble()
M_petsc.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)

# Zero out the rows corresponding to boundary elements
for i in range(0, DIMENSIONS):
    for j in range(0, MAT_LEN):
        A_petsc.setValue(top_elements[i], j, 0, addv=PETSc.InsertMode.INSERT)
        A_petsc.setValue(bottom_elements[i], j, 0, addv=PETSc.InsertMode.INSERT)
        M_petsc.setValue(top_elements[i], j, 0, addv=PETSc.InsertMode.INSERT)
        M_petsc.setValue(bottom_elements[i], j, 0, addv=PETSc.InsertMode.INSERT)

# Impose constant boundary conditions on the top and bottom boundaries (separately for each)
for i in range(0, DIMENSIONS):
    # Make the elements corresponding to a boundary element interacting with itself 1
    A_petsc.setValue(top_elements[i], top_elements[i], 1, addv=PETSc.InsertMode.INSERT)
    A_petsc.setValue(bottom_elements[i], bottom_elements[i], 1, addv=PETSc.InsertMode.INSERT)

    # Make the elements corresponding to a boundary element interacting with the one to its right -1
    if (i != DIMENSIONS-1):
        A_petsc.setValue(top_elements[i], top_elements[i+1], -1, addv=PETSc.InsertMode.INSERT)
        A_petsc.setValue(bottom_elements[i], bottom_elements[i+1], -1, addv=PETSc.InsertMode.INSERT)
    elif (i == DIMENSIONS-1):
        A_petsc.setValue(top_elements[DIMENSIONS-1], top_elements[0], -1, addv=PETSc.InsertMode.INSERT)
        A_petsc.setValue(bottom_elements[DIMENSIONS-1], bottom_elements[0], -1, addv=PETSc.InsertMode.INSERT)

# Re-assemble matrices to apply our changes
A_petsc.assemble()
M_petsc.assemble()

# Convert the matrices to a type that eigs can work with
A = sp.csr_matrix(A.mat().getValuesCSR()[::-1])
M = sp.csr_matrix(M.mat().getValuesCSR()[::-1])

# Get the first few eigenvalues and eigenvectors. Eigs will get the eigenvalues nearest to the value "sigma"
v, VV = eigs(A, 10, M, sigma=0)
print(v)
2

Complex eigenvalues are indicative of non-symmetric matrices. Are your final A and M symmetric? Based on your problem set-up, they should be. If they aren’t I’d try to figure out in which line of your code this symmetry is broken.

  • By the way, in what sense is this related to eigenvalues on the sphere? Your domain is a rectange, and in your formulation I am not recognizing transformations to spherical coordinates
  • You’re using PETSc for the assembly, but not SLEPc for the eigenvalues. Might be something to consider.