Mode solving for bent waveguides

Mode solving for bent waveguides#

reference values from [4]

Caution

This example is under construction eigenvalues match paper, overlaps not yet precise, working on it

As Julia-Dicts are not ordered, the mesh might become incorrect when adjusted (for now, better do the meshing in python)

using PyCall
np = pyimport("numpy")
shapely = pyimport("shapely")
shapely.affinity = pyimport("shapely.affinity")
clip_by_rect = pyimport("shapely.ops").clip_by_rect
OrderedDict = pyimport("collections").OrderedDict
mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict



function write_mesh(;
    radius = 2,
    wg_width = 0.5,
    wg_thickness = 0.22,
    sim_left = 0.5,
    sim_right = 4,
    sim_top = 1,
    sim_bottom = 4,
    pml_thickness = 3,
)
    core =
        shapely.geometry.box(radius - wg_width / 2, 0, radius + wg_width / 2, wg_thickness)

    env = shapely.geometry.box(
        radius - wg_width / 2 - sim_left,
        -sim_bottom - pml_thickness,
        radius + wg_width / 2 + sim_right + pml_thickness,
        wg_thickness + sim_top,
    )

    polygons = OrderedDict(
        core = core,
        box = clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
        clad = clip_by_rect(env, -np.inf, 0, np.inf, np.inf),
    )

    resolutions = Dict(
        "core" => Dict("resolution" => 0.01, "distance" => 1.5),
        "slab" => Dict("resolution" => 0.01, "distance" => 1.5),
    )

    mesh_from_OrderedDict(
        polygons,
        resolutions,
        default_resolution_max = 0.2,
        filename = "mesh.msh",
    )
end
Hide code cell output
write_mesh (generic function with 1 method)
using Gridap
using Gridap.Geometry
using Gridap.Visualization
using Gridap.ReferenceFEs
using GridapGmsh
using GridapMakie, CairoMakie

using Femwell.Maxwell.Waveguide

CairoMakie.inline!(true)
Hide code cell output
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_ckyULk".
 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/xG4T9/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/xG4T9/src/gl_backend.jl:2
 [18] include(mod::Module, _path::String)
    @ Base ./Base.jl:495
 [19] include(x::String)
    @ GLMakie ~/.julia/packages/GLMakie/xG4T9/src/GLMakie.jl:1
 [20] top-level scope
    @ ~/.julia/packages/GLMakie/xG4T9/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/xG4T9/src/gl_backend.jl:1
in expression starting at /home/runner/.julia/packages/GLMakie/xG4T9/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::IJulia.IJuliaStdio{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::IJulia.IJuliaStdio{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
radiuss = 1:0.5:5
wg_width = 0.5
sim_right = 4
sim_bottom = 4
pml_thickness = 3
neffs = ComplexF64[]
for radius in radiuss
    write_mesh(
        radius = radius,
        wg_width = wg_width,
        sim_right = sim_right,
        sim_bottom = sim_bottom,
    )
    model = GmshDiscreteModel("mesh.msh")
    Ω = Triangulation(model)
    labels = get_face_labeling(model)

    epsilons = ["core" => 3.48^2, "box" => 1.46^2, "clad" => 1.0^2]
    ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons)[tag]

    τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
    # pml_x = x -> 0.1 * max(0, x[1] - (radius + wg_width / 2 + sim_right))^2
    pml_order = 3
    pml_x =
        x ->
            1 / pml_order *
            max(0, (x[1] - (radius + wg_width / 2 + sim_right)) / pml_thickness)^pml_order *
            5
    pml_y = x -> 0.0

    modes = calculate_modes(model, ε  τ, λ = 1.55, order = 1)
    println(n_eff(modes[1]))
    modes = calculate_modes(
        model,
        ε  τ,
        λ = 1.55,
        order = 1,
        radius = radius,
        pml = [pml_x, pml_y],
        k0_guess = modes[1].k,
    )
    println(n_eff_cylidrical(modes[1]) / radius)
    println(log10(abs(imag(n_eff_cylidrical(modes[1]) / radius))))
    plot_mode(modes[1], absolute = false)
    push!(neffs, n_eff_cylidrical(modes[1]) / radius)
end

display(neffs)
radiuss_reference = 1:5
neff_fd = [2.40762, 2.39421, 2.39204, 2.39128, 2.39091]
log10imags_fd = [-3.68456, -6.41594, -9.37884, -9.98148, -10.39617]
neff_fmm = [2.40791, 2.39433, 2.39224, 2.39142, 2.39093]
log10imags_fmm = [-3.63721, -6.48982, -9.30488, -9.97048, -10.36615]
figure = Figure()
ax = Axis(figure[1, 1], xlabel = "Radius / μm", ylabel = "log(-imag(n_eff))")
lines!(ax, radiuss, log10.(-imag(neffs)))
scatter!(ax, radiuss, log10.(-imag(neffs)), label = "FEM")
scatter!(ax, radiuss_reference, log10imags_fd, label = "FD")
scatter!(
    ax,
    radiuss_reference,
    log10imags_fmm,
    label = "FMM",
    color = :red,
    marker = :xcross,
)
axislegend()
display(figure)
figure = Figure()
ax = Axis(figure[1, 1], xlabel = "Radius / μm", ylabel = "real(n_eff)")
lines!(ax, radiuss, real(neffs))
scatter!(ax, radiuss, real(neffs), label = "FEM")
scatter!(ax, radiuss_reference, neff_fd, label = "FD")
scatter!(ax, radiuss_reference, neff_fmm, label = "FMM", color = :red, marker = :xcross)
axislegend()
display(figure)