This tutorial is part of multiphenicsx.

Tutorial 04: first eigenvalue of the Laplacian with Dirichlet BCs

In this tutorial we evaluate the first eigenvalue of the Laplacian with homogeneous Dirichlet boundary conditions. Let $\Omega$ be the unit ball in 2D: given a constant $\alpha \in \mathbb{R}^+$, the goal is the find the smallest eigenvalue $\eta \in \mathbb{R}^+$ such that $$\begin{cases} -\alpha \Delta u = \eta u, & \text{in } \Omega,\\ u = 0, & \text{on } \partial\Omega, \end{cases}$$

The weak formulation of the eigenvalue problem is \begin{align*} &\text{find } (\eta, u) \in \mathbb{R} \times V \text{ s.t. }&\\ &\alpha \int_\Omega \nabla u \cdot \nabla v = \eta \int_\Omega \; u \; v, & \forall v \in V,\\ \end{align*} where $$ V = H^1_0(\Omega). $$

In the following, we will adopt the notation $\eta = \eta^{(\alpha)}$ when interested in comparing the first eigenvalue for different values of the parameter $\alpha$.

import typing

import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import scipy.special
import slepc4py.SLEPc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Geometrical parameters

mesh_size = 0.05

Mesh

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)
p1 = gmsh.model.geo.addPoint(0.0, 1.0, 0.0, mesh_size)
p2 = gmsh.model.geo.addPoint(0.0, -1.0, 0.0, mesh_size)
c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)
boundary = gmsh.model.geo.addCurveLoop([c0, c1])
domain = gmsh.model.geo.addPlaneSurface([boundary])

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [c0, c1], 1)
gmsh.model.addPhysicalGroup(2, [boundary], 0)
gmsh.model.mesh.generate(2)

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Circle)
Info    : [ 60%] Meshing curve 2 (Circle)
Info    : Done meshing 1D (Wall 0.00033851s, CPU 0.000659s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0338415s, CPU 0.033906s)
Info    : 1547 nodes 3093 elements

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

facets_partial_Omega = boundaries.indices[boundaries.values == 1]

viskex.dolfinx.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Approach #1

We define a $\mathbb{P}^2$ finite element space. In the simplest approach, we then use dolfinx to assemble the left-hand side matrix $A = A^{(\alpha)}$ representing the discretization of the Laplace operator, and the right-hand side matrix $B$ representing the discretization of the $L^2(\Omega)$ inner product. We then apply Dirichlet boundary conditions through dolfinx, and compare the obtained first eigenvalue with the analytical solution.

# Define a function space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))

# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Auxiliary function for the computation of the first eigenvalue using dolfinx
def get_smallest_eigenvalue_dolfinx(  # type: ignore[no-any-unimported]
    alpha: petsc4py.PETSc.RealType, diagonal_A: petsc4py.PETSc.RealType = 1.0
) -> tuple[
    petsc4py.PETSc.RealType, dolfinx.fem.Function
]:
    """Get the smallest eigenvalue, and the corresponding eigenfunction, using dolfinx."""
    # Define problem forms
    alpha_constant = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(alpha))
    a = alpha_constant * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    b = ufl.inner(u, v) * ufl.dx
    # Define boundary conditions
    zero = petsc4py.PETSc.ScalarType(0)
    bdofs_V = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, facets_partial_Omega)
    bc = dolfinx.fem.dirichletbc(zero, bdofs_V, V)
    # Assemble lhs and rhs matrices
    A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a), bcs=[bc], diagonal=diagonal_A)
    A.assemble()
    B = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(b), bcs=[bc])
    B.assemble()
    # Solve
    eps = slepc4py.SLEPc.EPS().create(mesh.comm)
    eps.setOperators(A, B)
    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GHEP)
    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)
    eps.getST().getKSP().setType("preonly")
    eps.getST().getKSP().getPC().setType("lu")
    eps.getST().getKSP().getPC().setFactorSolverType("mumps")
    eps.solve()
    assert eps.getConverged() >= 1
    # Extract first eigenvalue and eigenvector
    vr = dolfinx.la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
    vi = dolfinx.la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
    eigv = eps.getEigenpair(0, vr, vi)
    r, i = eigv.real, eigv.imag
    assert abs(i) < 1.e-10
    assert r > 0., "r = " + str(r) + " is not positive"
    # Destroy EPS object
    eps.destroy()
    # Transform eigenvector into an eigenfunction that can be plotted
    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    vr_fun = dolfinx.fem.Function(V)
    with vr_fun.x.petsc_vec.localForm() as vr_fun_local, \
            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, V.dofmap) as vr_wrapper:
        vr_fun_local[:] = vr_wrapper
    # Return eigenvalue and eigenfunction
    return r, vr_fun

In order to test this approach, recall that it can be shown by separation of variables in polar coordinates that $$\eta^{(1)} = j_{0,1}^2$$ where $j_{n, k}$ is the $k$-th positive zero of the $n$-th Bessel function $J_n$. We store the value of $j_{0,1}$ for later use.

j01,  = scipy.special.jn_zeros(0, 1)

j01, j01**2

(np.float64(2.4048255576957724), np.float64(5.783185962946783))

Once $\eta^{(1)}$ is known, the eigenvalue associated to any $\alpha$ can be easily obtained as $$ \eta^{(\alpha)} = \alpha j_{0,1}^2 $$ by linearity.

We first test this approach by computing the eigenvalue for the case $\alpha=0.1$.

eigenvalue_0p1, eigenfunction_0p1 = get_smallest_eigenvalue_dolfinx(0.1)

print(f"Computed: {eigenvalue_0p1}, vs expected {0.1 * j01**2}")

Computed: 0.5785616092123537, vs expected 0.5783185962946783

assert np.isclose(eigenvalue_0p1, 0.578561)

The approximation of $\eta^{(0.1)}$ is reasonably accurate. We can also plot the corresponding eigenfunction.

viskex.dolfinx.plot_scalar_field(eigenfunction_0p1, "eigenfunction 0.1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

We next try the case $\alpha = 1$.

eigenvalue_1, eigenfunction_1 = get_smallest_eigenvalue_dolfinx(1.0)

print(f"Computed: {eigenvalue_1}, vs expected {j01**2}")

Computed: 0.9999999999970641, vs expected 5.783185962946783

assert np.isclose(eigenvalue_1, 1.0)

The current approximation of $\eta^{(1)}$ is completely inaccurate. From the plot of the corresponding eigenfunction, we realize that the computed eigenvalue is a spurious one, which does not even correspond to an eigenfunction that satisfies the boundary conditions.

viskex.dolfinx.plot_scalar_field(eigenfunction_1, "eigenfunction 1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

boundary_integral_1 = mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(eigenfunction_1 * ufl.ds)), op=mpi4py.MPI.SUM)
print(f"Boundary integral of the eigenfunction: {boundary_integral_1}")
assert not np.isclose(boundary_integral_1, 0.)

Boundary integral of the eigenfunction: -0.33106000939833075

For comparison, the boundary integral of the eigenfunction associated to $\eta^{(0.1)}$ was indeed numerically zero.

boundary_integral_0p1 = mesh.comm.allreduce(
    dolfinx.fem.assemble_scalar(dolfinx.fem.form(eigenfunction_0p1 * ufl.ds)), op=mpi4py.MPI.SUM)
print(f"Boundary integral of the eigenfunction: {boundary_integral_0p1}")
assert np.isclose(boundary_integral_0p1, 0.)

Boundary integral of the eigenfunction: -2.377341490537274e-15

Approach #1: fixup

In the case $\alpha = 1$, the previous computation shows a spurious eigenvalue. The value resulting from there is a consequence of the way Dirichlet boundary conditions are internally applied by dolfinx in A = dolfinx.fem.petsc.assemble_matrix(..., bcs=[bc]): the row/columns of A associated to every DOF on $\partial\Omega$ are cleared up, and the diagonal is set to 1. By repeating the same procedure the matrix B on the right-hand side, we end up introducing a spurious eigenvalue equal to 1.

To obtain correct results, we should instead set the diagonal value of $A$ for row/columns associated to boundary DOFs to a number $d \in \mathbb{R}^+$, while still keeping a diagonal value equal to 1 for B. This will still introduce a spurious eigenvalue, with value $d$.

Let us try again for $\alpha = 1$, setting $d = 1.5$. We then expect the computed eigenvalue to be equal to $d$.

eigenvalue_1_fixup1, _ = get_smallest_eigenvalue_dolfinx(1.0, diagonal_A=1.5)

print(f"Computed: {eigenvalue_1_fixup1}, vs expected {j01**2}")

Computed: 1.5000000000037506, vs expected 5.783185962946783

assert np.isclose(eigenvalue_1_fixup1, 1.5)

We have successfully moved the first eigenvalue to $1.5$. Since we expect the true eigenvalue $\eta^{(1)}$ to be around 5.78, the goal will be to choose $d$ large enough, so that the spurious eigenvalue does not get returned as the minimum one. For istance, choose $d = 10$.

eigenvalue_1_fixup2, eigenfunction_1_fixup2 = get_smallest_eigenvalue_dolfinx(1.0, diagonal_A=10)

print(f"Computed: {eigenvalue_1_fixup2}, vs expected {j01**2}")

Computed: 5.785616092132607, vs expected 5.783185962946783

assert np.isclose(eigenvalue_1_fixup2, 5.785609)

A correct approximation of $\eta^{(1)}$ is now obtained. Furthermore, due to linear algebra properties we expect the first eigenfunction of the case $\alpha = 0.1$ and $\alpha = 1$ to be the same, up to a sign.

viskex.dolfinx.plot_scalar_field(eigenfunction_1_fixup2, "eigenfunction 1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

In order to compare the eigenfunctions associated to the two different values of $\alpha$, the next function normalizes eigenfunctions to ensure a consistent sign.

def normalize(u: dolfinx.fem.Function) -> None:
    """Normalize an eigenvector."""
    scaling_operations: list[tuple[  # type: ignore[no-any-unimported]
        dolfinx.fem.Function, typing.Callable[[dolfinx.fem.Function], ufl.Form],
        typing.Callable[[petsc4py.PETSc.ScalarType], petsc4py.PETSc.ScalarType]
    ]] = [
        # Scale functions with a W^{1,1} norm to take away possible sign differences.
        (u, lambda u: (u + u.dx(0) + u.dx(1)) * ufl.dx, lambda x: x),
        # Normalize functions with a H^1 norm.
        (u, lambda u: (ufl.inner(u, u) + ufl.inner(ufl.grad(u), ufl.grad(u))) * ufl.dx, lambda x: np.sqrt(x))
    ]
    for (function, bilinear_form, postprocess) in scaling_operations:
        scalar = postprocess(mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(bilinear_form(function))), op=mpi4py.MPI.SUM))
        function.x.petsc_vec.scale(1. / scalar)
        function.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)

normalize(eigenfunction_0p1)
normalize(eigenfunction_1_fixup2)

viskex.dolfinx.plot_scalar_field(eigenfunction_0p1, "eigenfunction 0.1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(eigenfunction_1_fixup2, "eigenfunction 1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

def compute_errors(u: dolfinx.fem.Function, uu: dolfinx.fem.Function) -> None:
    """Compute errors between two different cases."""
    u_norm_form = (ufl.inner(u, u) + ufl.inner(ufl.grad(u), ufl.grad(u))) * ufl.dx
    u_norm = np.sqrt(mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(u_norm_form)), op=mpi4py.MPI.SUM))
    err_norm_form = (ufl.inner(u - uu, u - uu) + ufl.inner(ufl.grad(u - uu), ufl.grad(u - uu))) * ufl.dx
    err_norm = np.sqrt(mesh.comm.allreduce(
        dolfinx.fem.assemble_scalar(dolfinx.fem.form(err_norm_form)), op=mpi4py.MPI.SUM))
    print("Relative error is equal to", err_norm / u_norm)
    assert np.isclose(err_norm / u_norm, 0., atol=1.e-6)

compute_errors(eigenfunction_0p1, eigenfunction_1_fixup2)

Relative error is equal to 5.69057624745014e-10

Approach #2

While the fix introduced above is attractive because of its simplicity, it has a potential drawback: in cases where no analytical solution exists, the user has to try different values of $d$ until they found a value which is suitable for shifting away the spurious eigenvalues.

In a second approach we leverage multiphenicsx capabilities to throw away DOFs associated to Dirichlet boundary conditions while assemblying A and B, so that they do not interfere with the eigenvalue calculation.

# Auxiliary function for the computation of the first eigenvalue using multiphenicsx
def get_smallest_eigenvalue_multiphenicsx(  # type: ignore[no-any-unimported]
    alpha: petsc4py.PETSc.RealType
) -> tuple[
    petsc4py.PETSc.RealType, dolfinx.fem.Function
]:
    """Get the smallest eigenvalue, and the corresponding eigenfunction, using multiphenicsx."""
    # Define restrictions.
    dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)
    dofs_V_partial_Omega = dolfinx.fem.locate_dofs_topological(V, boundaries.dim, facets_partial_Omega)
    restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, np.setdiff1d(dofs_V, dofs_V_partial_Omega))
    # Define problem forms
    alpha_constant = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(alpha))
    a = alpha_constant * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    b = ufl.inner(u, v) * ufl.dx
    # Assemble lhs and rhs matrices
    A = multiphenicsx.fem.petsc.assemble_matrix(
        dolfinx.fem.form(a), bcs=[], restriction=(restriction_V, restriction_V))
    A.assemble()
    B = multiphenicsx.fem.petsc.assemble_matrix(
        dolfinx.fem.form(b), bcs=[], restriction=(restriction_V, restriction_V))
    B.assemble()
    # Solve
    eps = slepc4py.SLEPc.EPS().create(mesh.comm)
    eps.setOperators(A, B)
    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GHEP)
    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)
    eps.getST().getKSP().setType("preonly")
    eps.getST().getKSP().getPC().setType("lu")
    eps.getST().getKSP().getPC().setFactorSolverType("mumps")
    eps.solve()
    assert eps.getConverged() >= 1
    # Extract first eigenvalue and eigenvector
    vr = dolfinx.la.create_petsc_vector(restriction_V.index_map, restriction_V.index_map_bs)
    vi = dolfinx.la.create_petsc_vector(restriction_V.index_map, restriction_V.index_map_bs)
    eigv = eps.getEigenpair(0, vr, vi)
    r, i = eigv.real, eigv.imag
    assert abs(i) < 1.e-10
    assert r > 0., "r = " + str(r) + " is not positive"
    # Destroy EPS object
    eps.destroy()
    # Transform eigenvector into an eigenfunction that can be plotted
    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    vr_fun = dolfinx.fem.Function(V)
    with vr_fun.x.petsc_vec.localForm() as vr_fun_local, \
            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, V.dofmap, restriction_V) as vr_wrapper:
        vr_fun_local[:] = vr_wrapper
    # Return eigenvalue and eigenfunction
    return r, vr_fun

We test this approach on both cases $\alpha = 0.1$ and $\alpha = 1$.

eigenvalue_0p1_approach2, eigenfunction_0p1_approach2 = get_smallest_eigenvalue_multiphenicsx(0.1)

print(f"Computed: {eigenvalue_0p1}, vs expected {0.1 * j01**2}")

Computed: 0.5785616092123537, vs expected 0.5783185962946783

assert np.isclose(eigenvalue_0p1, 0.578561)

eigenvalue_1_approach2, eigenfunction_1_approach2 = get_smallest_eigenvalue_multiphenicsx(1.0)

print(f"Computed: {eigenvalue_1_approach2}, vs expected {j01**2}")

Computed: 5.785616092128737, vs expected 5.783185962946783

assert np.isclose(eigenvalue_1_approach2, 5.785609)

We see that the computed values match the expected ones for $\eta^{(0.1)}$ and $\eta^{(1)}$. We also compare eigenfunctions between the two approaches.

normalize(eigenfunction_0p1_approach2)
normalize(eigenfunction_1_approach2)

viskex.dolfinx.plot_scalar_field(eigenfunction_0p1_approach2, "eigenfunction 0.1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(eigenfunction_1_approach2, "eigenfunction 1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

compute_errors(eigenfunction_0p1, eigenfunction_0p1_approach2)

Relative error is equal to 8.479722274403712e-07

compute_errors(eigenfunction_0p1, eigenfunction_1_approach2)

Relative error is equal to 8.479151066453871e-07