Helmholtz eigenvalue problem in a rectangle

Helmholtz eigenvalue problem in a rectangle#

\[\begin{split} \mathbb{S} \begin{cases} \Omega = [0, L_x] \times [0, L_y] \\ u_{\text{D}}\vert_{\partial\Omega}=0 \\ f = 0 \\ \end{cases} \end{split}\]
import numpy as np
from dolfinx.fem import FunctionSpace

from lucifex.mesh import rectangle_mesh, mesh_boundary
from lucifex.solver import (
    evp, BoundaryConditions, OptionsSLEPc,
)
from lucifex.viz import plot_colormap, save_figure
from lucifex.pde.eigen import helmholtz


Lx = 1.0
Ly = 1.0
Nx = 32
Ny = 32

mesh = rectangle_mesh(Lx, Ly, Nx, Ny)
boundary = mesh_boundary(
        mesh, 
        {
            "left": lambda x: x[0],
            "right": lambda x: x[0] - Lx,
            "lower": lambda x: x[1],
            "upper": lambda x: x[1] - Ly,
        },
)
bcs = BoundaryConditions(('dirichlet', boundary.union, 0.0))
fs = FunctionSpace(mesh, ('P', 1))

nev = 10
slepc = OptionsSLEPc(
    eps_tol=1e-10,
    eps_target=0.0,
    eps_nev=nev,
    eps_ncv=50,
    eps_max_it=1000,
    eps_which='smallest_real',
)
u_solver = evp(helmholtz, bcs, slepc)(fs)
u_solver.solve()
sqrt = lambda x: np.sqrt(np.real(x)) if np.real(x) > 0 else None
k_values = [sqrt(-e) for e in u_solver.eigenvalues]
k_values
[None,
 None,
 None,
 4.444667393076721,
 7.034411761557017,
 7.034411761557114,
 8.900046596911208,
 9.967342410168534,
 9.967342410168753,
 11.361152904419761]
for k, u in zip(k_values, u_solver.eigenfunctions):
    if k is not None:
        fig, ax = plot_colormap(u, title=f'$k={k}$')
        save_figure(f'k={k:.3f}\nu', thumbnail=(k is k_values[-1]))(fig)
../../_images/91ed16fb9b0cad976ef48d7207386aa11afc11b599f1d28567b3dde60315cff1.png ../../_images/0b96976d0e39dbb21dbdac73d7050242657c422f068db959cdd089943523f83a.png ../../_images/86abd646ccef14cc7c3674538b6b160f63bf344c7df59a2d4ec3d812e583cbac.png ../../_images/f6c28c3728e0fb08ed3b7f60ada755e462cedfb9a0715683ee5f51736777b3a5.png ../../_images/0b49830526acfa54686002f06de9bd0e1fa19d1f996b4424c4726ba7f0300891.png ../../_images/2c320a2cd29720b7f0af1a9c13f05d84b127d3a61d02c75342a0a61d7db65914.png ../../_images/86be75a2dc24907d0518a95158b71462622f59a1936ca35a2afb4fe74b71ea8b.png
The Kernel crashed while executing code in the current cell or a previous cell. 

Please review the code in the cell(s) to identify a possible cause of the failure. 

Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. 

View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.