diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst index 1f4bf8f6c2..fb6c299d1e 100644 --- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst +++ b/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 `__. -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 `__ matrices -and eigenvalue solvers in `SLEPc `__. +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. diff --git a/docs/source/conf.py b/docs/source/conf.py index 47f140d76e..8670269803 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -262,6 +262,7 @@ \sphinxDUC{2218}{$\circ$} \sphinxDUC{22C5}{$\cdot$} \sphinxDUC{25A3}{$\boxdot$} +\sphinxDUC{03BB}{$\lambda$} % Sphinx equivalent of % \DeclareUnicodeCharacter{}{} diff --git a/docs/source/team.ini b/docs/source/team.ini index 12388d8be9..f0b79413db 100644 --- a/docs/source/team.ini +++ b/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: diff --git a/firedrake/__init__.py b/firedrake/__init__.py index 13a1a46f49..369fb85274 100644 --- a/firedrake/__init__.py +++ b/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 * diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py new file mode 100644 index 0000000000..c883bf96be --- /dev/null +++ b/firedrake/eigensolver.py @@ -0,0 +1,194 @@ +"""Specify and solve finite element eigenproblems.""" +from firedrake.assemble import assemble +from firedrake.function import Function +from firedrake import utils +from firedrake.petsc import OptionsManager, flatten_parameters +from firedrake.exceptions import ConvergenceError +try: + from slepc4py import SLEPc +except ImportError: + SLEPc = None +__all__ = ["LinearEigenproblem", + "LinearEigensolver"] + + +class LinearEigenproblem(): + """Generalised linear eigenvalue problem. + + The problem has the form, find `u`, `λ` such that:: + + A(u, v) = λM(u, v) ∀ v ∈ V + + Parameters + ---------- + A : ufl.Form + The bilinear form A(u, v). + M : ufl.Form + The mass form M(u, v), defaults to inner(u, v) * dx. + bcs : DirichletBC or list of DirichletBC + The boundary conditions. + bc_shift: float + The value to shift the boundary condition eigenvalues by. + + Notes + ----- + + If Dirichlet boundary conditions are supplied then these will result in the + eigenproblem having a nullspace spanned by the basis functions with support + on the boundary. To facilitate solution, this is shifted by the specified + amount. It is the user's responsibility to ensure that the shift is not + close to an actual eigenvalue of the system. + """ + def __init__(self, A, M=None, bcs=None, bc_shift=0.0): + if not SLEPc: + raise ImportError( + "Unable to import SLEPc, eigenvalue computation not possible " + "(try firedrake-update --slepc)" + ) + + self.A = A # LHS + args = A.arguments() + v, u = args + if M: + self.M = M + else: + from ufl import inner, dx + self.M = inner(u, v) * dx + self.output_space = u.function_space() + self.bcs = bcs + self.bc_shift = bc_shift + + def dirichlet_bcs(self): + """Return an iterator over the Dirichlet boundary conditions.""" + for bc in self.bcs: + yield from bc.dirichlet_bcs() + + @utils.cached_property + def dm(self): + r"""Return the dm associated with the output space.""" + return self.output_space.dm + + +class LinearEigensolver(OptionsManager): + r"""Solve a LinearEigenproblem. + + Parameters + ---------- + problem : LinearEigenproblem + The eigenproblem to solve. + n_evals : int + The number of eigenvalues to compute. + options_prefix : str + The options prefix to use for the eigensolver. + solver_parameters : dict + PETSc options for the eigenvalue problem. + ncv : int + Maximum dimension of the subspace to be used by the solver. See + `SLEPc.EPS.setDimensions`. + mpd : int + Maximum dimension allowed for the projected problem. See + `SLEPc.EPS.setDimensions`. + + Notes + ----- + + Users will typically wish to set solver parameters specifying the symmetry + of the eigenproblem and which eigenvalues to search for first. + + The former is set using the options available for `EPSSetProblemType + `__. + + For example if the bilinear form is symmetric (Hermitian in complex mode), + one would add this entry to `solver_options`:: + + "eps_gen_hermitian": None + + As always when specifying PETSc options, `None` indicates that the option + in question is a flag and hence doesn't take an argument. + + The eigenvalues to search for first are specified using the options for + `EPSSetWhichEigenPairs `__. + + For example, to look for the eigenvalues with largest real part, one + would add this entry to `solver_options`:: + + "eps_largest_real": None + """ + + DEFAULT_EPS_PARAMETERS = {"eps_type": "krylovschur", + "eps_tol": 1e-10, + "eps_target": 0.0} + + def __init__(self, problem, n_evals, *, options_prefix=None, + solver_parameters=None, ncv=None, mpd=None): + + self.es = SLEPc.EPS().create(comm=problem.dm.comm) + self._problem = problem + self.n_evals = n_evals + self.ncv = ncv + self.mpd = mpd + solver_parameters = flatten_parameters(solver_parameters or {}) + for key in self.DEFAULT_EPS_PARAMETERS: + value = self.DEFAULT_EPS_PARAMETERS[key] + solver_parameters.setdefault(key, value) + if self._problem.bcs: + solver_parameters.setdefault("st_type", "sinvert") + super().__init__(solver_parameters, options_prefix) + self.set_from_options(self.es) + + def check_es_convergence(self): + r"""Check the convergence of the eigenvalue problem.""" + r = self.es.getConvergedReason() + try: + reason = SLEPc.EPS.ConvergedReasons[r] + except KeyError: + reason = ("unknown reason (petsc4py enum incomplete?), " + "try with -eps_converged_reason") + if r < 0: + raise ConvergenceError( + r"""Eigenproblem failed to converge after %d iterations. + Reason: + %s""" % (self.es.getIterationNumber(), reason) + ) + + def solve(self): + r"""Solve the eigenproblem. + + Returns + ------- + int + The number of Eigenvalues found. + """ + self.A_mat = assemble(self._problem.A, bcs=self._problem.bcs).M.handle + self.M_mat = assemble( + self._problem.M, bcs=self._problem.bcs, + weight=self._problem.bc_shift and 1./self._problem.bc_shift + ).M.handle + + self.es.setDimensions(nev=self.n_evals, ncv=self.ncv, mpd=self.mpd) + self.es.setOperators(self.A_mat, self.M_mat) + with self.inserted_options(): + self.es.solve() + nconv = self.es.getConverged() + if nconv == 0: + raise ConvergenceError("Did not converge any eigenvalues.") + return nconv + + def eigenvalue(self, i): + r"""Return the i-th eigenvalue of the solved problem.""" + return self.es.getEigenvalue(i) + + def eigenfunction(self, i): + r"""Return the i-th eigenfunction of the solved problem. + + Returns + ------- + (Function, Function) + The real and imaginary parts of the eigenfunction. + """ + eigenmodes_real = Function(self._problem.output_space) # fn of V + eigenmodes_imag = Function(self._problem.output_space) + with eigenmodes_real.dat.vec_wo as vr: + with eigenmodes_imag.dat.vec_wo as vi: + self.es.getEigenvector(i, vr, vi) # gets the i-th eigenvector + return eigenmodes_real, eigenmodes_imag # returns Firedrake fns diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py new file mode 100644 index 0000000000..f58339238f --- /dev/null +++ b/tests/regression/test_eigensolver.py @@ -0,0 +1,82 @@ +import numpy as np +import pytest +from firedrake import * + + +def evals(n, degree=1, mesh=None): + '''We base this test on the 1D Poisson problem with Dirichlet boundary + conditions, outlined in part 1 of Daniele Boffi's + 'Finite element approximation of eigenvalue problems' Acta Numerica 2010''' + # Create mesh and define function space + if mesh is None: + mesh = IntervalMesh(n, 0, pi) + V = FunctionSpace(mesh, "CG", degree) + # Define variational problem + u = TrialFunction(V) + v = TestFunction(V) + a = (inner(grad(u), grad(v))) * dx + + # Create eigenproblem with boundary conditions + bc = DirichletBC(V, 0.0, "on_boundary") + eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-6666.) + + # Create corresponding eigensolver, looking for n eigenvalues + eigensolver = LinearEigensolver( + eigenprob, n, solver_parameters={"eps_largest_real": None} + ) + ncov = eigensolver.solve() + + # boffi solns + h = pi / n + true_values = np.zeros(ncov-1) + estimates = np.zeros(ncov-1) + for k in range(ncov-1): + true_val = 6 / h**2 + # k+1 because we skip the trivial 0 eigenvalue + true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h)) + true_values[k] = true_val + + estimates[k] = eigensolver.eigenvalue(k).real + + true_values[-1] = eigenprob.bc_shift + + # sort in case order of numerical and analytic values differs. + return sorted(true_values), sorted(estimates) + + +@pytest.mark.parametrize(('n', 'degree', 'tolerance'), + [(5, 1, 1e-13), + (10, 1, 1e-13), + (20, 1, 1e-13), + (30, 1, 1e-13)]) +def test_evals_1d(n, degree, tolerance): + true_values, estimates = evals(n, degree=degree) + assert np.allclose(true_values, estimates, rtol=tolerance) + + +def poisson_eigenvalue_2d(i): + mesh = RectangleMesh(10*2**i, 10*2**i, pi, pi) + V = FunctionSpace(mesh, "CG", 1) + u = TrialFunction(V) + v = TestFunction(V) + + bc = DirichletBC(V, 0.0, "on_boundary") + + ep = LinearEigenproblem(inner(grad(u), grad(v)) * dx, + bcs=bc, bc_shift=666.0) + + es = LinearEigensolver(ep, 1, solver_parameters={"eps_gen_hermitian": None, + "eps_largest_real": None}) + + es.solve() + return es.eigenvalue(0)-2.0 + + +def test_evals_2d(): + """2D Eigenvalue convergence test. As with Boffi, we observe that the + convergence rate convergest to 2 from above.""" + errors = np.array([poisson_eigenvalue_2d(i) for i in range(5)]) + + convergence = np.log(errors[:-1]/errors[1:])/np.log(2.0) + + assert all(convergence > 2.0)