Eigensolver (#2985)

Introduce LinearEigenproblem and LinearEigensolver classes. Also add tests and a demo.

---------

Co-authored-by: David Ham <[email protected]>
Co-authored-by: ksagiyam <[email protected]>

Author: 3 people

demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst

@@ -6,20 +6,16 @@ Oceanic Basin Modes: Quasi-Geostrophic approach
    This tutorial was contributed by Christine Kaufhold and `Francis
    Poulin <mailto:[email protected]>`__.
 
-As a continuation of the Quasi-Geostrophic (QG) model described in the
-other tutorial, we will now see how we can use Firedrake to compute
-the spatial structure and frequencies of the freely evolving modes in this system,
-what are referred to as basin modes.
-Oceanic basin modes are low frequency structures that propagate
-zonally in the oceans that alter the dynamics of Western Boundary Currents,
-such as the Gulf Stream. In this particular tutorial we will show how to
-solve the QG eigenvalue problem with no basic state and no dissipative
-forces.
-Unlike the other demo that integrated the equations forward in time, in this
-problem it is necessary to compute the eigenvalues and eigenfunctions
-for a particular differential operator. This requires using
-`PETSc <http://www.mcs.anl.gov/petsc/>`__ matrices
-and eigenvalue solvers in `SLEPc <http://slepc.upv.es>`__.
+As a continuation of the Quasi-Geostrophic (QG) model described in the other
+tutorial, we will now see how we can use Firedrake to compute the spatial
+structure and frequencies of the freely evolving modes in this system, what are
+referred to as basin modes. Oceanic basin modes are low frequency structures
+that propagate zonally in the oceans that alter the dynamics of Western
+Boundary Currents, such as the Gulf Stream. In this particular tutorial we will
+show how to solve the QG eigenvalue problem with no basic state and no
+dissipative forces. Unlike the other demo that integrated the equations forward
+in time, in this problem it is necessary to compute the eigenvalues and
+eigenfunctions for a particular differential operator.
 
 This demo requires SLEPc and slepc4py to be installed.  This is most easily
 achieved by providing the optional `--slepc` flag to either `firedrake-install`
@@ -105,133 +101,97 @@ Firedrake code
 --------------
 
 Using this form, we can now implement this eigenvalue problem in
-Firedrake. We import the Firedrake, PETSc, and SLEPc libraries. ::
-
-   from firedrake import *
-   from firedrake.petsc import PETSc
-   try:
-       from slepc4py import SLEPc
-   except ImportError:
-       import sys
-       warning("Unable to import SLEPc, eigenvalue computation not possible (try firedrake-update --slepc)")
-       sys.exit(0)
+Firedrake. We start by importing Firedrake. ::
 
+  from firedrake import *
 
 We specify the geometry to be a square geometry with :math:`50` cells
 with length :math:`1`. ::
 
-   Lx   = 1.
-   Ly   = 1.
-   n0   = 50
-   mesh = RectangleMesh(n0, n0, Lx, Ly, reorder=None)
+  Lx   = 1.
+  Ly   = 1.
+  n0   = 50
+  mesh = RectangleMesh(n0, n0, Lx, Ly, reorder=None)
 
 Next we define the function spaces within which our solution will
 reside. ::
 
-   Vcg  = FunctionSpace(mesh,'CG',3)
+  Vcg  = FunctionSpace(mesh,'CG',3)
 
 We impose zero Dirichlet boundary conditions, in a strong sense, which
 guarantee that we have no-normal flow at the boundary walls. ::
 
-   bc = DirichletBC(Vcg, 0.0, "on_boundary")
+  bc = DirichletBC(Vcg, 0.0, "on_boundary")
 
 The two non-dimensional parameters are the :math:`\beta` parameter, set
 by the sphericity of the Earth, and the Froude number, the relative
 importance of rotation to stratification. ::
 
-   beta = Constant('1.0')
-   F    = Constant('1.0')
-
-Additionally, we can create some Functions to store the eigenmodes. ::
-
-   eigenmodes_real, eigenmodes_imag = Function(Vcg), Function(Vcg)
+  beta = Constant('1.0')
+  F    = Constant('1.0')
 
 We define the Test Function :math:`\phi` and the Trial Function
 :math:`\psi` in our function space. ::
 
-   phi, psi = TestFunction(Vcg), TrialFunction(Vcg)
+  phi, psi = TestFunction(Vcg), TrialFunction(Vcg)
 
 To build the weak formulation of our equation we need to build two PETSc
 matrices in the form of a generalized eigenvalue problem,
-:math:`A\psi = \lambda M\psi`. We impose the boundary conditions on the
-mass matrix :math:`M`, since that is where we used integration by parts. ::
-
-   a =  beta*phi*psi.dx(0)*dx
-   m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx
-   petsc_a = assemble(a).M.handle
-   petsc_m = assemble(m, bcs=bc).M.handle
+:math:`A\psi = \lambda M\psi` ::
 
-We can declare how many eigenpairs, eigenfunctions and eigenvalues, we
-want to find ::
+  eigenproblem = LinearEigenproblem(
+          A=beta*phi*psi.dx(0)*dx,
+          M=-inner(grad(psi), grad(phi))*dx - F*psi*phi*dx,
+          bcs=bc)
 
-   num_eigenvalues = 1
+Next we program our eigenvalue solver through the PETSc options system. The
+first is specifying that we have an generalized eigenvalue problem that is
+nonhermitian. Then, we ask for the eigenvalues with the largest imaginary
+part. Finally we set the spectral transform to shift with no target::
 
-Next we will impose parameters onto our eigenvalue solver. The first is
-specifying that we have an generalized eigenvalue problem that is
-nonhermitian. The second specifies the spectral transform shift factor
-to be non-zero. The third requires we use a Krylov-Schur method,
-which is the default so this is not strictly necessary. Then, we ask for
-the eigenvalues with the largest imaginary part. Finally, we specify the
-tolerance. ::
+  opts = {"eps_gen_non_hermitian": None,
+          "eps_largest_imaginary": None,
+          "st_type": "shift",
+          "eps_target": None,
+          "st_pc_factor_shift_type": "NONZERO"}
 
-   opts = PETSc.Options()
-   opts.setValue("eps_gen_non_hermitian", None)
-   opts.setValue("st_pc_factor_shift_type", "NONZERO")
-   opts.setValue("eps_type", "krylovschur")
-   opts.setValue("eps_largest_imaginary", None)
-   opts.setValue("eps_tol", 1e-10)
+Finally, we build our eigenvalue solver, specifying in this case that we just
+want to see the first eigenvalue, eigenvector pair::
 
-Finally, we build our eigenvalue solver using SLEPc. We add our PETSc
-matrices into the solver as operators and use setFromOptions() to call
-the PETSc parameters we previously declared. ::
+  eigensolver = LinearEigensolver(eigenproblem, n_evals=1,
+                                  solver_parameters=opts)
 
-   es = SLEPc.EPS().create(comm=COMM_WORLD)
-   es.setDimensions(num_eigenvalues)
-   es.setOperators(petsc_a, petsc_m)
-   es.setFromOptions()
-   es.solve()
+Now solve the system. This returns the number of converged eigenvalues. ::
 
-Additionally we can find the number of converged eigenvalues. ::
-
-   nconv = es.getConverged()
+  nconv = eigensolver.solve()
 
 We now get the real and imaginary parts of the eigenvalue and
 eigenvector for the leading eigenpair (that with the largest in
-magnitude imaginary part).  First we check if we actually managed to
-converge any eigenvalues at all. ::
-
-   if nconv == 0:
-       import sys
-       warning("Did not converge any eigenvalues")
-       sys.exit(0)
-
-If we did, we go ahead and extract them from the SLEPc eigenvalue
-solver::
-
-   vr, vi = petsc_a.getVecs()
+magnitude imaginary part). ::
 
-   lam = es.getEigenpair(0, vr, vi)
+  lam = eigensolver.eigenvalue(0)
 
-and we gather the final eigenfunctions ::
+and we gather the corresponding eigenfunctions ::
 
-   eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi
+  eigenmode_real, eigenmode_imag = eigensolver.eigenfunction(0)
 
 We can now list and show plots for the eigenvalues and eigenfunctions
 that were found. ::
 
-   print("Leading eigenvalue is:", lam)
+  print("Leading eigenvalue is:", lam)
 
-   try:
-       import matplotlib.pyplot as plt
-       fig, axes = plt.subplots()
-       colors = tripcolor(eigenmodes_real, axes=axes)
-       fig.colorbar(colors)
+  try:
+      import matplotlib.pyplot as plt
+      fig, axes = plt.subplots()
+      colors = tripcolor(eigenmode_real, axes=axes)
+      fig.colorbar(colors)
 
-       fig, axes = plt.subplots()
-       colors = tripcolor(eigenmodes_imag, axes=axes)
-       fig.colorbar(colors)
-   except ImportError:
-       warning("Matplotlib not available, not plotting eigemodes")
+      fig, axes = plt.subplots()
+      colors = tripcolor(eigenmode_imag, axes=axes)
+      fig.colorbar(colors)
+      plt.show()
+  except ImportError:
+      warning("Matplotlib not available, not plotting eigemodes")
 
 Below is a plot of the spatial structure of the real part of one of the eigenmodes computed above.
 

docs/source/conf.py

@@ -262,6 +262,7 @@
 \sphinxDUC{2218}{$\circ$}
 \sphinxDUC{22C5}{$\cdot$}
 \sphinxDUC{25A3}{$\boxdot$}
+\sphinxDUC{03BB}{$\lambda$}
 % Sphinx equivalent of
 % \DeclareUnicodeCharacter{}{}
 

docs/source/team.ini

@@ -111,4 +111,5 @@ Ben Sepanski: https://bensepanski.github.io
 Jemma Shipton:
 Joseph G. Wallwork: https://www.imperial.ac.uk/people/j.wallwork16
 Florian Wechsung: https://florianwechsung.github.io
+Yian Zeng:
 Fangyi Zhou:

firedrake/__init__.py

@@ -98,6 +98,7 @@
 from firedrake.ufl_expr import *
 from firedrake.utility_meshes import *
 from firedrake.variational_solver import *
+from firedrake.eigensolver import *
 from firedrake.vector import *
 from firedrake.version import __version__ as ver, __version_info__, check  # noqa: F401
 from firedrake.ensemble import *