Benchmark of the mode solver 2

Reproducing [1], where the modes of a strip and several rib waveguide were calculated and presented with an error value. The error for all modes is calculated to be smaller than \(\pm 3 \cdot 10^{-6}\), thus this should be the maximum derivation for our simulations. The results are presented here:

Show code cell source Hide code cell source

from collections import OrderedDict

import numpy as np
import pandas as pd
import shapely
import shapely.affinity
from juliacall import Main as jl
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

julia_script = [
    "using Gridap",
    "using Gridap.Geometry",
    "using GridapGmsh",
    "using Femwell.Maxwell.Waveguide",
    "using GridapMakie, CairoMakie",
    "import PhysicalConstants.CODATA2018: c_0, μ_0, ε_0",
    "import Unitful: ustrip",
]

for line in julia_script:
    jl.seval(line)

core_width = 3
core_thickness = 1
slab_width = 18
box_thickness = 10
clad_thickness = 3

slab_thicknesses = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7]
neff_values_paper = [
    3.412022,
    3.412126,
    3.412279,
    3.412492,
    3.412774,
    3.413132,
    3.413571,
    3.414100,
]
neff_values_femwell_slepc = []
neff_values_femwell_scipy = []
neff_values_julia = []

for slab_thickness in slab_thicknesses:
    if slab_thickness == 0.0:
        core = shapely.Polygon(
            (
                (-core_width / 2, 0),
                (-core_width / 2, core_thickness),
                (core_width / 2, core_thickness),
                (core_width / 2, 0),
            )
        )
    else:
        core = shapely.Polygon(
            (
                (-slab_width / 2, 0),
                (-slab_width / 2, slab_thickness),
                (-core_width / 2, slab_thickness),
                (-core_width / 2, core_thickness),
                (core_width / 2, core_thickness),
                (core_width / 2, slab_thickness),
                (slab_width / 2, slab_thickness),
                (slab_width / 2, 0),
            )
        )

    polygons = OrderedDict(
        core=core,
        box=shapely.box(-slab_width / 2, -box_thickness, slab_width / 2, 0),
        clad=shapely.box(-slab_width / 2, clad_thickness, slab_width / 2, 0),
    )

    resolutions = dict(core={"resolution": 0.1, "distance": 3})

    mesh = from_meshio(
        mesh_from_OrderedDict(
            polygons, resolutions, default_resolution_max=0.3, filename="mesh.msh"
        )
    )
    refractive_indices = {"core": 3.44, "box": 3.40, "clad": 1}
    jl.refractive_indices = refractive_indices

    julia_script = [
        'model = GmshDiscreteModel("mesh.msh")',
        "Ω = Triangulation(model)",
        "labels = get_face_labeling(model)",
        "τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)",
        'refractive_indices = Dict("core" => 3.44^2, "box" => 3.40^2, "clad" => 1^2)',
        "ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in refractive_indices)[tag]",
        "modes = calculate_modes(model, ε ∘ τ,  λ = 1.15, order = 1)",
        "neffs = [real(n_eff(mode)) for mode in modes]",
    ]
    for line in julia_script:
        jl.seval(line)

    neff_values_julia.append(jl.neffs[0])

    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros()
    for subdomain, n in refractive_indices.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2

    # modes = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2, solver="slepc")
    # neff_values_femwell_slepc.append(np.real(modes[0].n_eff))

    modes = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2, solver="scipy")
    neff_values_femwell_scipy.append(np.real(modes[0].n_eff))

Show code cell output Hide code cell output

[juliapkg] Found dependencies: /home/runner/miniconda3/lib/python3.12/site-packages/juliacall/juliapkg.json
[juliapkg] Found dependencies: /home/runner/miniconda3/lib/python3.12/site-packages/juliapkg/juliapkg.json
[juliapkg] Locating Julia ~1.6.1, ~1.7, ~1.8, ~1.9, =1.10.0, ^1.10.3
[juliapkg] Using Julia 1.10.3 at /opt/hostedtoolcache/julia/1.10.3/x64/bin/julia
[juliapkg] Using Julia project at /home/runner/work/femwell/femwell

Detected IPython. Loading juliacall extension. See https://juliapy.github.io/PythonCall.jl/stable/compat/#IPython

---------------------------------------------------------------------------
JuliaError                                Traceback (most recent call last)
Cell In[1], line 25
     14 julia_script = [
     15     "using Gridap",
     16     "using Gridap.Geometry",
   (...)
     21     "import Unitful: ustrip",
     22 ]
     24 for line in julia_script:
---> 25     jl.seval(line)
     27 core_width = 3
     28 core_thickness = 1

File ~/.julia/packages/PythonCall/Nr75f/src/JlWrap/module.jl:27, in seval(self, expr)
     25         return ValueBase.__dir__(self) + self._jl_callmethod($(pyjl_methodnum(pyjlmodule_dir)))
     26     def seval(self, expr):
---> 27         return self._jl_callmethod($(pyjl_methodnum(pyjlmodule_seval)), expr)
     28 """,
     29             @__FILE__(),

JuliaError: The following 1 direct dependency failed to precompile:

GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a]

Failed to precompile GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a] to "/home/runner/.julia/compiled/v1.10/GLMakie/jl_YTvOvl".
┌ Warning:     OpenGL/GLFW wasn't loaded correctly or couldn't be initialized.
│     This likely means, you're on a headless server without having OpenGL support setup correctly.
│     Have a look at the troubleshooting section in the readme:
│     https://github.com/MakieOrg/Makie.jl/tree/master/GLMakie#troubleshooting-opengl.
└ @ GLMakie ~/.julia/packages/GLMakie/o5Nad/src/gl_backend.jl:4
ERROR: LoadError: InitError: Exception[GLFW.GLFWError(65550, "X11: The DISPLAY environment variable is missing"), ErrorException("glfwInit failed")]
Stacktrace:
  [1] __init__()
    @ GLFW ~/.julia/packages/GLFW/wmoTL/src/GLFW.jl:35
  [2] run_module_init(mod::Module, i::Int64)
    @ Base ./loading.jl:1134
  [3] register_restored_modules(sv::Core.SimpleVector, pkg::Base.PkgId, path::String)
    @ Base ./loading.jl:1122
  [4] _include_from_serialized(pkg::Base.PkgId, path::String, ocachepath::String, depmods::Vector{Any})
    @ Base ./loading.jl:1067
  [5] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128)
    @ Base ./loading.jl:1581
  [6] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:1938
  [7] __require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1812
  [8] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
  [9] invoke_in_world
    @ ./essentials.jl:923 [inlined]
 [10] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1803
 [11] macro expansion
    @ ./loading.jl:1790 [inlined]
 [12] macro expansion
    @ ./lock.jl:267 [inlined]
 [13] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1753
 [14] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
 [15] invoke_in_world
    @ ./essentials.jl:923 [inlined]
 [16] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1746
 [17] top-level scope
    @ ~/.julia/packages/GLMakie/o5Nad/src/gl_backend.jl:2
 [18] include(mod::Module, _path::String)
    @ Base ./Base.jl:495
 [19] include(x::String)
    @ GLMakie ~/.julia/packages/GLMakie/o5Nad/src/GLMakie.jl:1
 [20] top-level scope
    @ ~/.julia/packages/GLMakie/o5Nad/src/GLMakie.jl:87
 [21] include
    @ ./Base.jl:495 [inlined]
 [22] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::Nothing)
    @ Base ./loading.jl:2222
 [23] top-level scope
    @ stdin:3
during initialization of module GLFW
in expression starting at /home/runner/.julia/packages/GLMakie/o5Nad/src/gl_backend.jl:1
in expression starting at /home/runner/.julia/packages/GLMakie/o5Nad/src/GLMakie.jl:1
in expression starting at stdin:
Stacktrace:
  [1] pkgerror(msg::String)
    @ Pkg.Types /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/Types.jl:70
  [2] precompile(ctx::Pkg.Types.Context, pkgs::Vector{Pkg.Types.PackageSpec}; internal_call::Bool, strict::Bool, warn_loaded::Bool, already_instantiated::Bool, timing::Bool, _from_loading::Bool, kwargs::@Kwargs{io::Base.PipeEndpoint})
    @ Pkg.API /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:1659
  [3] precompile(pkgs::Vector{Pkg.Types.PackageSpec}; io::Base.PipeEndpoint, kwargs::@Kwargs{_from_loading::Bool})
    @ Pkg.API /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:159
  [4] precompile
    @ /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:147 [inlined]
  [5] #precompile#114
    @ /opt/hostedtoolcache/julia/1.10.3/x64/share/julia/stdlib/v1.10/Pkg/src/API.jl:146 [inlined]
  [6] #invokelatest#2
    @ ./essentials.jl:894 [inlined]
  [7] invokelatest
    @ ./essentials.jl:889 [inlined]
  [8] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:1963
  [9] __require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1812
 [10] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
 [11] invoke_in_world
    @ ./essentials.jl:923 [inlined]
 [12] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1803
 [13] macro expansion
    @ ./loading.jl:1790 [inlined]
 [14] macro expansion
    @ ./lock.jl:267 [inlined]
 [15] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1753
 [16] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
 [17] invoke_in_world
    @ ./essentials.jl:923 [inlined]
 [18] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1746
 [19] eval
    @ ./boot.jl:385 [inlined]
 [20] eval
    @ ./Base.jl:88 [inlined]
 [21] pyjlmodule_seval(self::Module, expr::Py)
    @ PythonCall.JlWrap ~/.julia/packages/PythonCall/Nr75f/src/JlWrap/module.jl:13
 [22] _pyjl_callmethod(f::Any, self_::Ptr{PythonCall.C.PyObject}, args_::Ptr{PythonCall.C.PyObject}, nargs::Int64)
    @ PythonCall.JlWrap ~/.julia/packages/PythonCall/Nr75f/src/JlWrap/base.jl:67
 [23] _pyjl_callmethod(o::Ptr{PythonCall.C.PyObject}, args::Ptr{PythonCall.C.PyObject})
    @ PythonCall.JlWrap.Cjl ~/.julia/packages/PythonCall/Nr75f/src/JlWrap/C.jl:63

Show code cell source Hide code cell source

pd.DataFrame(
    {
        "slab_thickness": slab_thicknesses,
        "reference value": (f"{n:.6f}" for n in neff_values_paper),
        # "calculated value slepc": (f"{n:.6f}" for n in neff_values_femwell_slepc),
        # "difference slepc": (
        #    f"{n1-n2:.6f}" for n1, n2 in zip(neff_values_paper, neff_values_femwell_slepc)
        # ),
        "calculated value scipy": (f"{n:.6f}" for n in neff_values_femwell_scipy),
        "difference scipy": (
            f"{n1-n2:.6f}" for n1, n2 in zip(neff_values_paper, neff_values_femwell_scipy)
        ),
        "calculated value julia": (f"{n:.6f}" for n in neff_values_julia),
        "difference julia": (
            f"{n1-n2:.6f}" for n1, n2 in zip(neff_values_paper, neff_values_julia)
        ),
    }
).style.apply(
    lambda differences: [
        "background-color: green" if abs(float(difference)) < 5e-6 else "background-color: red"
        for difference in differences
    ],
    subset=["difference scipy", "difference julia"],  # "difference slepc",
)

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.