Benchmark of the mode solver 1

Contents

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:

Hide 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:595, in compute_modes(basis_epsilon_r, epsilon_r, wavelength, mu_r, num_modes, order, metallic_boundaries, radius, n_guess, solver)
    592     sigma = sigma = k0**2 * np.max(epsilon_r) * 1.1
    594 if metallic_boundaries:
--> 595     lams, xs = solve(
    596         *condense(
    597             -A,
    598             -B,
    599             D=basis.get_dofs(None if metallic_boundaries is True else metallic_boundaries),
    600             x=basis.zeros(dtype=complex),
    601         ),
    602         solver=solver(k=num_modes, sigma=sigma),
    603     )
    604 else:
    605     lams, xs = solve(
    606         -A,
    607         -B,
    608         solver=solver(k=num_modes, sigma=sigma),
    609     )

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.

Contents