DOLFIN

eigenvalue problems with Dirichlet boundary conditions

Asked by Doug Arnold

This concerns the computation of eigenvalues, e.g., for the Laplacian with Dirichlet boundary conditions. This is tricky, since in dolfin we work with matrices whose size is the number of DOFs including those on the boundary, but number of eigenvalues for the eigenvalue problem we are considering is only the number of DOFs excluding those on the Dirichlet boundary. Looking at dolfin email, I see that people are unclear about (1) how to do this mathematically and (2) how to implement in it dolfin. I have figured out the solution to (1), and the implemention in the case we run serially, but don't see how to implement it implement it in parallel.

This is important since it is natural to expect dolfin to solve Dirichlet eigenvalue problems in parallel. If anyone has done this and could share some code please do. Meanwhile, I'll say what I have done. At the end is a question, which I hope someone can help me with.

(1) Suppose that A is the stiffness matrix and M is the mass matrix, so the eigenvalue problem with Neumann boundary conditions is the (A,M) generalized eigenvalue problem

  A u = lambda M u,

which can be solved in FEniCS using SLEPc.

With Dirichlet boundary conditions on some or all of the boundary, the problem should be changed by

 (a) zeroing out the rows and columns of A corresponding to boundary nodes, but putting a 1 (or other nonzero) on the diagonal of such rows;
 (b) zeroing out the rows and columns of M corresponding to boundary node including the diagonal elements.

With these changes the eigenpairs to the new (A,M) generalized eigenvalue problem are exactly the eigenpairs for the problem with Dirichlet boundary conditions. (Because M is singular, the number of eigenvalues is less than the size of the matrices.)

(2) Now to implement this in dolfin, we can make the matrix A as follows:

   a = dot(grad(u), grad(v))*dx
   b = v*dx
   bc = DirichletBC(V, Constant(0.0), bdry)
   A, B = assemble_system(a, b, bc)

(Remarks: (1) b and B are of no further use. (2) we have to set A = down_cast(A) to convert A to a PETSc matrix before calling SLEPc.)

The problem comes with constructing the matrix M in dolfin. I don't know a way to simply set the rows and columns associated to boundary DOFs to zero without putting something on the diagonal. So what I do now is compute the list of indices of the boundary DOFs and then use

  for i in bdrydoflist:
    row = np.array((i,), dtype='I')
    col = row
    val = np.array( (0.,) )
    M.set(val, row, col)
    M.apply("insert")

This seems ugly, but it works fine in serial. Not surprisingly, it does not work in parallel.

All this leads to my question: is there a nice way to construct the matrix M described in (b) above. Or, is there way to set the diagonal elements corresponding to boundary DOFs to 0? PETSc has MatDiagonalSet for this, but I don't know how to use this in dolfin.

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Joachim Haga
Solved:
Last query:
Last reply:
Revision history for this message
Best Joachim Haga (jobh) said : 		#1

> All this leads to my question: is there a nice way to construct the matrix M described in (b) above. Or, is there way to set the diagonal elements corresponding to boundary DOFs to 0? PETSc has MatDiagonalSet for this, but I don't know how to use this in dolfin.

I think it should work to construct a second matrix containing just
the boundary, and subtract it from the original. Something along these
lines:

A, B = assemble_system(a, b, bcs)
C, D = assemble_system(Constant(0)*a, b, bcs)

A0 = A-C

(
You can also use symmetric_assemble instead of assemble_system,
A, A_asymm = symmetric_assemble(a, bcs)
...
)

-j.

Revision history for this message
Doug Arnold (dnarnold) said : 		#2

That's a great idea, Joachim, and it works fine.

If anyone wants to see a simple code to solve the Dirichlet eigenvalue problem let me know.

>> You can also use symmetric_assemble instead of assemble_system,

I've been looking for something like this, but I guess it was introduced after version 1.0.0.
I don't find it in my version.

Thanks!

Revision history for this message
Joachim Haga (jobh) said : 		#3

>>> You can also use symmetric_assemble instead of assemble_system,
>
> I've been looking for something like this, but I guess it was introduced after version 1.0.0.
> I don't find it in my version.

That's true, it's a later addition.

I just remembered that there is a function bc.zero(A) that is a more
direct approach. It is normally not symmetry preserving, but when the
rows are diagonal (when Dirichlet BCs are already applied
symmetrically to the rows) it is, of course.

-j.

Revision history for this message
Doug Arnold (dnarnold) said : 		#4

That's even better. So all one needs to do to impose Dirichlet bcs on an eigenvalue problem is

  A = symmetric_assemble(a, bc)
  M = symmetric_assemble(m, bc)
  bc.zero(M)

and it works in parallel too.

Revision history for this message
Joachim Haga (jobh) said : 		#5

> Doug Arnold posted a new comment:
> That's even better. So all one needs to do to impose Dirichlet bcs on
> an eigenvalue problem is
>
> A = symmetric_assemble(a, bc)
> M = symmetric_assemble(m, bc)
> bc.zero(M)
>
> and it works in parallel too.

Yes. Almost. However:

- symmetric_assemble is not available in 1.0 (as you've seen), and I
just realised it's disabled/scheduled for removal even there. Will fix
the latter problem; for the former, assemble_system can be used
instead.

- symmetric_assemble always returns symmetric AND antisymmetric
matrices, you'll have to throw away the antisymmetric part:

A, _ = symmetric_assemble(a, bc)
M, _ = symmetric_assemble(m, bc)
bc.zero(M)

- I haven't worked on this in a while, you'll need to check that
zero() actually does what I think it does.

Revision history for this message
Doug Arnold (dnarnold) said : 		#6

I had already checked that it works with assemble_system. Here is a complete eigenvalue code
which gives the right answer in serial and parallel. Thanks a lot for your help.

# This program demonstrates the use of SLEPc to
# find the three smallest magnitude eigenvalues of the Laplacian
# with Dirichlet boundary conditions on a square and the corresponding
# eigenfunctions.

from dolfin import *
import numpy as np

# Test for PETSc and SLEPc
if not has_linear_algebra_backend("PETSc"):
    print "DOLFIN has not been configured with PETSc. Exiting."
    exit()

if not has_slepc():
    print "DOLFIN has not been configured with SLEPc. Exiting."
    exit()

def geteig(n, deg, nreq, export_eigenfunctions):
  """
  Compute the eigenvalues of the Laplacian with Dirichlet boundary conditions on the square.
  Use a mesh of n x n squares, divided into triangles, Lagrange elements of degree deg, and
  request nreq eigenpairs. If export_eigenfunctions=True, write the eigenfunctions to
  PVD files. Return values are the number of converged eigenpairs and the computed eigenvalues.
  """

  # Define basis and bilinear form
  mesh = UnitSquare(n, n)
  V = FunctionSpace(mesh, "CG", deg)
  u = TrialFunction(V)
  v = TestFunction(V)
  a = dot(grad(u), grad(v))*dx
  m = u*v*dx

  # define Dirichlet boundary conditions
  def bdry(x, on_boundary):
    return on_boundary

  bc = DirichletBC(V, Constant(0.0), bdry)

  # Assemble the stiffness matrix and the mass matrix.
  b = v*dx # just need this to feed an argument to assemble_system
  A, _ = assemble_system(a, b, bc)
  M, _ = assemble_system(m, b, bc)
  # set the diagonal elements of M corresponding to boundary nodes to zero to
  # remove spurious eigenvalues.
  bc.zero(M)
  # downcast to PETSc matrices
  A = down_cast(A)
  M = down_cast(M)

  # Create eigensolver
  eigensolver = SLEPcEigenSolver(A, M)
  # Specify the solution method (default is krylov-schur)
  eigensolver.parameters["solver"] = "krylov-schur"
  # Specify the part of the spectrum desired
  eigensolver.parameters["spectrum"] = "smallest magnitude"
  # Specify the problem type (this can make a big difference)
  eigensolver.parameters["problem_type"] = "gen_hermitian"
  # Use the shift-and-invert spectral transformation
  eigensolver.parameters["spectral_transform"] = "shift-and-invert"
  # Specify the shift
  eigensolver.parameters["spectral_shift"] = 1.0e-10

  # Compute the eigenvalues. Here is where we specify the number of eigenvalues.
  nreq = 3
  eigensolver.solve(nreq)

  # Check the number of eigenvalues that converged.
  # Extract the eigenvalues (ignore the imaginary part) and compare with the exact values
  nconv = eigensolver.get_number_converged()
  if MPI.process_number() == 0:
    print "Number of eigenvalues successfully computed: ", nconv
  eigenvalues = []
  for i in range(nconv):
    r, c, rx, cx = eigensolver.get_eigenpair(i) # real and complex part of eigenvalue
    eigenvalues.append(r)
    # export the eigenfunctions for external visualization
    if export_eigenfunctions:
      eigenfun = Function(V)
      eigenfun.vector()[:] = rx
      File("eig{0}.pvd".format(i)) << eigenfun
  return nconv, eigenvalues

# main program
exacteigenvalues = (2*pi*pi, 5*pi*pi, 5*pi*pi)
n = 25; deg = 2; nreq = 3; export_eigenfunctions = True
nconv, eigenvalues = geteig(n, deg, nreq, export_eigenfunctions)
if MPI.process_number() == 0:
  print "Smallest positive eigenvalues computed and exact: "
  for i in range(min(nconv, nreq)):
    print eigenvalues[i], exacteigenvalues[i]

Revision history for this message
Garth Wells (garth-wells) said : 		#7

On Tue, Oct 2, 2012 at 3:35 PM, Joachim Haga
<email address hidden> wrote:
> Question #210106 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210106
>
> Joachim Haga posted a new comment:
>> Doug Arnold posted a new comment:
>> That's even better. So all one needs to do to impose Dirichlet bcs on
>> an eigenvalue problem is
>>
>> A = symmetric_assemble(a, bc)
>> M = symmetric_assemble(m, bc)
>> bc.zero(M)
>>
>> and it works in parallel too.
>
> Yes. Almost. However:
>
> - symmetric_assemble is not available in 1.0 (as you've seen), and I
> just realised it's disabled/scheduled for removal even there. Will fix
> the latter problem; for the former, assemble_system can be used
> instead.
>
> - symmetric_assemble always returns symmetric AND antisymmetric
> matrices, you'll have to throw away the antisymmetric part:
>
> A, _ = symmetric_assemble(a, bc)
> M, _ = symmetric_assemble(m, bc)
> bc.zero(M)
>

The latest system_assembler will (shortly) be able to:

   A0 = assemble(a)
   A1 = assemble(a, bcs)

Garth

> - I haven't worked on this in a while, you'll need to check that
> zero() actually does what I think it does.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Revision history for this message
Doug Arnold (dnarnold) said : 		#8

> The latest system_assembler will (shortly) be able to:
>
> A0 = assemble(a)
> A1 = assemble(a, bcs)

It would be nice if it were possible to control show the bcs are applied to A1. Sometimes
you want to zero out the rows corresponding to boundary DOFs, sometimes the rows
and columns, and in either case, sometimes you want to put a nonzero on the diagonals.
It would be nice if there were a clear way to do that. For the eigenvalue
problem you want to zero out rows and colums, put nonzeros on the diagonal of the
stiffness matrix, but leave zeros on the diagonal of the mass matrix.