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)
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.