Benchmark of the mode solver 1#
Reproducing [1], where the modes of a analytically solvable geometry are calculated. The error for all modes is calculated to be smaller than \(\pm 1 \cdot 10^{-8}\). We’ll show that we get pretty close, but will stop at a resonable resolution to keep the runtime sensible. Getting even higher accurancy will be left open for adaptive refinement. The results are presented here:
Show code cell source
from collections import OrderedDict
import numpy as np
import pandas as pd
import shapely
import shapely.affinity
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio
from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
epsilons_paper = [
{"core": 2.25, "clad": 1},
{"core": 8, "clad": 1},
{"core": 1, "clad": 2.25},
{"core": 1, "clad": 8},
]
boundaries = [["left", "right"], ["left", "right"], ["left", "top"], ["left", "top"]]
neff_values_paper = [1.27627404, 2.65679692, 1.387926425, 2.761465320]
neff_values_femwell_slepc = []
neff_values_femwell_scipy = []
polygons = OrderedDict(
left=shapely.LineString(((0, 0), (0, 1))),
right=shapely.LineString(((1, 0), (1, 1))),
top=shapely.LineString(((0, 1), (1, 1))),
core=shapely.box(0, 0, 0.5, 0.5),
clad=shapely.box(0, 0, 1, 1),
)
mesh = from_meshio(
mesh_from_OrderedDict(polygons, {}, default_resolution_max=0.01, filename="mesh.msh")
)
basis0 = Basis(mesh, ElementTriP0())
for epsilons, boundaries in zip(epsilons_paper, boundaries):
epsilon = basis0.zeros()
for subdomain, e in epsilons.items():
epsilon[basis0.get_dofs(elements=subdomain)] = e
modes = compute_modes(
basis0,
epsilon,
wavelength=1.5,
num_modes=1,
order=2,
metallic_boundaries=boundaries,
solver="slepc",
)
neff_values_femwell_slepc.append(np.real(modes[0].n_eff))
modes = compute_modes(
basis0,
epsilon,
wavelength=1.5,
num_modes=1,
order=2,
metallic_boundaries=boundaries,
solver="scipy",
)
neff_values_femwell_scipy.append(np.real(modes[0].n_eff))
pd.DataFrame(
{
"Epsilons": [
f"{epsilons['core']:.2f} / {epsilons['clad']:.2f}" for epsilons in epsilons_paper
],
"Reference value": (f"{n:.8f}" for n in neff_values_paper),
"Calculated value slepc": (f"{n:.8f}" for n in neff_values_femwell_slepc),
"Difference slepc": (
f"{n1-n2:.8f}" for n1, n2 in zip(neff_values_paper, neff_values_femwell_slepc)
),
"Calculated value scipy": (f"{n:.8f}" for n in neff_values_femwell_scipy),
"Difference scipy": (
f"{n1-n2:.8f}" for n1, n2 in zip(neff_values_paper, neff_values_femwell_scipy)
),
}
).style.apply(
lambda differences: [
"background-color: green" if abs(float(difference)) < 5e-6 else "background-color: red"
for difference in differences
],
subset=["Difference slepc", "Difference scipy"],
)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[1], line 43
40 for subdomain, e in epsilons.items():
41 epsilon[basis0.get_dofs(elements=subdomain)] = e
---> 43 modes = compute_modes(
44 basis0,
45 epsilon,
46 wavelength=1.5,
47 num_modes=1,
48 order=2,
49 metallic_boundaries=boundaries,
50 solver="slepc",
51 )
52 neff_values_femwell_slepc.append(np.real(modes[0].n_eff))
54 modes = compute_modes(
55 basis0,
56 epsilon,
(...) 61 solver="scipy",
62 )
File ~/miniconda3/lib/python3.12/site-packages/femwell/maxwell/waveguide.py:564, in compute_modes(basis_epsilon_r, epsilon_r, wavelength, mu_r, num_modes, order, metallic_boundaries, radius, n_guess, solver)
561 sigma = sigma = k0**2 * np.max(epsilon_r) * 1.1
563 if metallic_boundaries:
--> 564 lams, xs = solve(
565 *condense(
566 -A,
567 -B,
568 D=basis.get_dofs(None if metallic_boundaries is True else metallic_boundaries),
569 x=basis.zeros(dtype=complex),
570 ),
571 solver=solver(k=num_modes, sigma=sigma),
572 )
573 else:
574 lams, xs = solve(
575 -A,
576 -B,
577 solver=solver(k=num_modes, sigma=sigma),
578 )
File ~/miniconda3/lib/python3.12/site-packages/skfem/utils.py:268, in solve(A, b, x, I, solver, **kwargs)
266 logger.info("Solving linear system, shape={}.".format(A.shape))
267 if isinstance(b, spmatrix):
--> 268 out = solve_eigen(A, b, x, I, solver, **kwargs) # type: ignore
269 elif isinstance(b, ndarray):
270 out = solve_linear(A, b, x, I, solver, **kwargs) # type: ignore
File ~/miniconda3/lib/python3.12/site-packages/skfem/utils.py:211, in solve_eigen(A, M, x, I, solver, **kwargs)
208 solver = solver_eigen_scipy(**kwargs)
210 if x is not None and I is not None:
--> 211 L, X = solver(A, M, **kwargs)
212 y = np.tile(x.copy()[:, None], (1, X.shape[1]))
213 if isinstance(I, tuple):
File ~/miniconda3/lib/python3.12/site-packages/femwell/solver.py:93, in solver_eigen_slepc.<locals>.solver(K, M, **solve_time_kwargs)
92 def solver(K, M, **solve_time_kwargs):
---> 93 from petsc4py import PETSc
94 from slepc4py import SLEPc
96 which = {
97 "LM": SLEPc.EPS.Which.LARGEST_MAGNITUDE,
98 "SM": SLEPc.EPS.Which.SMALLEST_MAGNITUDE,
(...) 106 "ALL": SLEPc.EPS.Which.ALL,
107 }
File ~/miniconda3/lib/python3.12/site-packages/petsc4py/PETSc.py:4
1 ARCH = None
2 from petsc4py.lib import ImportPETSc # noqa: E402
----> 4 PETSc = ImportPETSc(ARCH)
5 PETSc._initialize()
6 del PETSc
File ~/miniconda3/lib/python3.12/site-packages/petsc4py/lib/__init__.py:30, in ImportPETSc(arch)
26 """
27 Import the PETSc extension module for a given configuration name.
28 """
29 path, arch = getPathArchPETSc(arch)
---> 30 return Import('petsc4py', 'PETSc', path, arch)
File ~/miniconda3/lib/python3.12/site-packages/petsc4py/lib/__init__.py:97, in Import(pkg, name, path, arch)
95 warnings.filterwarnings('ignore', message='numpy.ndarray size changed')
96 # import extension module from 'path/arch' directory
---> 97 module = import_module(pkg, name, path, arch)
98 module.__arch__ = arch # save arch value
99 setattr(sys.modules[pkg], name, module)
File ~/miniconda3/lib/python3.12/site-packages/petsc4py/lib/__init__.py:76, in Import.<locals>.import_module(pkg, name, path, arch)
74 module = importlib.util.module_from_spec(spec)
75 sys.modules[fullname] = module
---> 76 spec.loader.exec_module(module)
77 return module
78 f, fn, info = imp.find_module(name, pathlist)
File petsc4py/PETSc.pyx:1, in init petsc4py.PETSc()
ValueError: numpy.dtype size changed, may indicate binary incompatibility. Expected 96 from C header, got 88 from PyObject
Bibliography#
[1]
G.R. Hadley. High-accuracy finite-difference equations for dielectric waveguide analysis II: dielectric corners. Journal of Lightwave Technology, 20(7):1219–1231, July 2002. URL: https://doi.org/10.1109/jlt.2002.800371, doi:10.1109/jlt.2002.800371.