.. |fenics| replace:: FEniCS .. |slepc| replace:: ``SLEPc`` .. _slepc : https://slepc.upv.es/documentation/ .. _helmholtz: Helmholtz equation ================== This tutorial demonstrates how to solve the `Helmholtz equation `_ (the eigenvalue problem for the Laplace operator) on a box mesh with an opposite inlet and outlet. Specifically, it shows how to: * obtain the variational formulation of an eigenvalue problem * apply Dirichlet boundary conditions (tricky!) * use |slepc| and ``scipy`` to solve the generalized eigenvalue problem Additionally, different |slepc|_ solvers are compared and a benchmark is performed against AVSP (our in-house tool to solve this kind of problems). Formulating the problem ----------------------- Our goal is to find :math:`u \neq 0, \lambda \in \mathbb{R}` such that .. math:: \begin{aligned} &-\Delta u=\lambda u \text { in } \Omega \\ &u=0 \text { on } \partial \Omega_{D} \\ &\frac{\partial u}{\partial n} =0 \text { on } \partial \Omega_N \\ \end{aligned} Following the usual approach (remember FEniCS needs the variational formulation of a problem), we multiply both sides by a test function :math:`v` and integrate by parts. Our problem is to find :math:`0 \neq u \in V, \lambda \in \mathbb{R}` such that .. math:: \int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x=\lambda \int_{\Omega} u v \mathrm{~d} x \quad \forall v \in \hat{V} Notice our eigenvalue problem is now a generalized eigenvalue problem of the form: .. math:: A x=\lambda M x This is very different from the problems most commonly solved in FEniCS, as :math:`u` appears in both sides of the equation, and therefore the problem cannot be represented in the canonical notation for variational problems: find :math:`u \in V` such that .. math:: a(u, v)=L(v) \quad \forall v \in \hat{V} You can find more details about the formulation in `these `_ delightful notes from Patrick Farrell. A simple solution ----------------- Creating the mesh ***************** For this example we will create a simple box mesh. .. code-block:: python from dolfin.cpp.generation import BoxMesh from dolfin.cpp.generation import Point N = 60 lx, ly, lz = 0.2, 0.1, 1.0 mesh = BoxMesh(Point(0, 0, 0), Point(lx, ly, lz), int(lx*N), int(ly*N), int(lz*N)) This is how it looks: .. image:: images/mesh.png :align: center :scale: 80% Translate the problem formulation into code ******************************************* The equations shown above are reflected in the code as follows: .. code-block:: python from dolfin.function.functionspace import FunctionSpace from dolfin.function.argument import TrialFunction from dolfin.function.argument import TestFunction from ufl import inner from ufl import grad from ufl import dx V = FunctionSpace(mesh, 'Lagrange', 1) u_h = TrialFunction(V) v_h = TestFunction(V) a = inner(grad(u_h), grad(v_h)) * dx b = inner(u_h, v_h) * dx Assemble the matrices and apply boundary conditions *************************************************** The assembling of the matrices and application of the boundary conditions is the most sensitive part of this tutorial. The main reason is that `it is not advised `_ to remove rows and columns of matrices, as it is an inefficient operation (remember we are dealing with sparse matrices). To circumvent that, `several alternatives `_ to impose Dirichlet boundary conditions exist. Independently of the chosen strategy, the key point is to preserve symmetry in both matrices, as the available eigensolvers are better suited for symmetric problems (and it is a pitty to loose such a nice property!). Our strategy consists in zeroing-out the rows and columns corresponding to "Dirichlet degrees of freedom" and assign :math:`\alpha` and 1 to the corresponding diagonals, for :math:`A` and :math:`M`, respectively. :math:`\alpha` is the value that meaningless (artificial) eigenvalues will take. If you want to avoid them, just choose a value that is outside the desired spectrum. Let's now perform the assembly: .. code-block:: python from dolfin.cpp.la import PETScMatrix from dolfin.fem.assembling import assemble A = PETScMatrix(V.mesh().mpi_comm()) assemble(a, tensor=A) B = PETScMatrix(V.mesh().mpi_comm()) assemble(b, tensor=B) Notice we are using ``PETSc`` matrices, which are sparse. These data type is required to work with |slepc| later. For the application of boundary conditions, we start by collecting the Dirichlet dofs: .. code-block:: python bc_dofs = [] for bc in bcs: bc_dofs.extend(list(bc.get_boundary_values().keys())) Here, ``bcs`` is a list with Dirichlet boundary conditions. We are fixing the two faces perpendicular to :math:`z` axis. Finally, we use directly a method available in the ``PETSc`` matrix object to apply our chosen strategy (|fenics| equivalent method does not work in parallel): .. code-block:: python A.mat().zeroRowsColumnsLocal(bc_dofs, diag=diag_value) B.mat().zeroRowsColumnsLocal(bc_dofs, diag=1.) Note that homogeneous Neumann boundary conditions are very simple to treat: do nothing! What to do with these matrices? ******************************* |fenics| considers |slepc| the go-to library to solve eigenvalue problems. Nevertheless, after assembling the matrices, we can choose any solver we want (including our own!). A good starting point is the most commonly used scientific libraries ``numpy`` and ``scipy``. ``numpy`` provides machinery to solve "normal" eigenvalue problems, but it is not able to solve generalized eigenvalue problems. Besides, the lack of sparse data structures clearly shows its inability to scale (you cannot go that far with dense matrices). On the other hand, ``scipy`` provide solutions both for `dense `_ and `sparse `_ matrices. Let's solve it with ``scipy`` then: .. code-block:: import scipy A_array = scipy.sparse.csr_matrix(A.mat().getValuesCSR()[::-1]) B_array = scipy.sparse.csr_matrix(B.mat().getValuesCSR()[::-1]) k = 20 # number of eigenvalues which = 'SM' # smallest magnityde w, v = scipy.linalg.eigh(A_array, B_array, k=k, which=which) `scipy.sparse.linalg.eigsh `_ relies on `ARPACK's `_ implementation of the Implicitly Restarted Lanczos method. On the other hand, the basic usage of |slepc| is also very simple: .. code-block:: python import numpy as np from dolfin.cpp.la import SLEPcEigenSolver # instantiate solver solver = SLEPcEigenSolver(A, B) # update parameters solver.parameters['solver'] = 'krylov-schur' solver.parameters['spectrum'] = 'smallest magnitude' solver.parameters['problem_type'] = 'gen_hermitian' solver.parameters['tolerance'] = 1e-4 # solve n_eig = 20 solver.solve(n_eig) # collect eigenpairs w, v = [], [] for i in range(solver.get_number_converged()): r, _, rv, _ = solver.get_eigenpair(i) w.append(r) v.append(rv) w = np.array(w) v = np.array(v).T Our preliminary studies show ``scipy`` is not competitive with ``SLEPc`` for large problems. This `paper `_ presents open-source libraries capable of solving eigenvalue problems (focus is generalized non-hermitian eigenvalue problems), which may also be explored: LAPACK, ARPACK, Anasazi, SLEPc, FEAST, z-PARES. Their summary of each library (and algorithms) is worth reading. What `SLEPc` offers? ******************** Above, we've used |slepc| without much thinking. Yet, it offers several `algorithms `_ to solve eigenvalue problems, including: power iterations, subspace iterations, Arnoldi, Lanczos, Krylov-Schur (default), Jacobi-Davidson and Generalized Davidson. |fenics| documentation is outdated (it is a good idea to peek into the code to know what is really `available `_). These are the possible arguments for ``solver.parameters['solver']``: ``power``, ``subspace``, ``arnoldi``, ``lanczos``, ``krylov-schur``, ``lapack``, ``arpack``, ``jacobi-davidson``, ``generalized-davidson``. You can find `here `_ further information about each option. We are simply relying on the methods |fenics| exposes from |slepc| library. If we want to go further, we can also use `slepc4py `_, which are the bindings for |slepc| and provide much more flexibility. The following picture shows how each algorithm compare in terms of performance for the problem described above (changing the mesh size): .. image:: images/algo_comparison.png :align: center | As you can see, ``generalized-davidson`` outperforms all the other (the options not tested are not suitable for this problem). In fact, it allows to solve a problem with approximately 100k degrees of freedom in about 13 seconds (very dependent on parameters such as tolerance). From this number of degrees of freedom on, the weakest point of the chosen strategy is not anymore the solver, but the assembling step (solution: go parallel). It is also important to check the number of converged solutions (we ask for 20): .. image:: images/converged_comp.png :align: center Modes ***** And an example of the modes 1, 2, 6 and 7 visualized with ``paraview``: .. image:: images/mode1.png :align: center :scale: 40% | .. image:: images/mode2.png :align: center :scale: 40% | .. image:: images/mode6.png :align: center :scale: 40% | .. image:: images/mode7.png :align: center :scale: 40% | Closing remarks --------------- Have you seen how simple it is to solve an eigenproblem with |fenics|? An example script can be found `here `_. .. toctree:: :maxdepth: 1 :hidden: benchmark_avsp kokiy_fenics