Mode solving

Contents

Mode solving#

Caution

This example is under construction 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

wg_width = 0.5
wg_thickness = 0.22
slab_width = 3
slab_thickness = 0.11
core = shapely.geometry.box(-wg_width / 2, 0, wg_width / 2, wg_thickness)
slab = shapely.geometry.box(slab_width / 2, 0, slab_width / 2, slab_thickness)
env = shapely.affinity.scale(core.buffer(3, resolution = 8), xfact = 0.5)

polygons = OrderedDict(
    core = core,
    slab = slab,
    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.03, "distance" => 0.5),
    "slab" => Dict("resolution" => 0.03, "distance" => 0.5),
)

mesh = mesh_from_OrderedDict(
    polygons,
    resolutions,
    default_resolution_max = 1,
    filename = "mesh.msh",
)

Hide code cell output

 
PyObject <meshio mesh object>
  Number of points: 632
  Number of cells:
    vertex: 1
    vertex: 1
    line: 4
    line: 17
    line: 8
    line: 17
    line: 8
    line: 6
    line: 8
    line: 3
    triangle: 344
    triangle: 352
    triangle: 528
  Cell sets: slab_points, slab, core___box, core___clad, box___clad, core, box, clad, gmsh:bounding_entities
  Point data: gmsh:dim_tags
  Cell data: gmsh:physical, gmsh:geometrical
  Field data: slab_points, slab, core___box, core___clad, box___clad, core, box, clad

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

 
InitError: could not load library "/home/runner/.julia/artifacts/b2e490a224990d92e99cc4463b8df48c8c454971/lib/libpangocairo-1.0.so"
/home/runner/.julia/artifacts/b2e490a224990d92e99cc4463b8df48c8c454971/lib/libpangoft2-1.0.so.0: undefined symbol: FcConfigSetDefaultSubstitute
during initialization of module Pango_jll

Stacktrace:
  [1] dlopen(s::String, flags::UInt32; throw_error::Bool)
    @ Base.Libc.Libdl ./libdl.jl:117
  [2] dlopen(s::String, flags::UInt32)
    @ Base.Libc.Libdl ./libdl.jl:116
  [3] macro expansion
    @ ~/.julia/packages/JLLWrappers/m2Pjh/src/products/library_generators.jl:63 [inlined]
  [4] __init__()
    @ Pango_jll ~/.julia/packages/Pango_jll/CWpUN/src/wrappers/x86_64-linux-gnu.jl:22
  [5] run_module_init(mod::Module, i::Int64)
    @ Base ./loading.jl:1134
  [6] register_restored_modules(sv::Core.SimpleVector, pkg::Base.PkgId, path::String)
    @ Base ./loading.jl:1122
  [7] _include_from_serialized(pkg::Base.PkgId, path::String, ocachepath::String, depmods::Vector{Any})
    @ Base ./loading.jl:1067
  [8] _tryrequire_from_serialized(modkey::Base.PkgId, path::String, ocachepath::String, sourcepath::String, depmods::Vector{…})
    @ Base ./loading.jl:1481
  [9] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128)
    @ Base ./loading.jl:1574
 [10] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:1938
 [11] __require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1812
 [12] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
 [13] invoke_in_world
    @ ./essentials.jl:923 [inlined]
 [14] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1803
 [15] macro expansion
    @ ./loading.jl:1790 [inlined]
 [16] macro expansion
    @ ./lock.jl:267 [inlined]
 [17] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1753
 [18] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
 [19] invoke_in_world
    @ ./essentials.jl:923 [inlined]
 [20] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1746

model = GmshDiscreteModel("mesh.msh")
Ω = Triangulation(model)
#fig = plot(Ω)
#fig.axis.aspect=DataAspect()
#wireframe!(Ω, color=:black, linewidth=1)
#display(fig)

labels = get_face_labeling(model)

epsilons = ["core" => 3.5^2, "slab" => 1.44^2, "box" => 1.444^2, "clad" => 1.44^2]
ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons)[tag]


#dΩ = Measure(Ω, 1)
τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)

modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 1, order = 1)
println(n_eff(modes[1]))
# write_mode_to_vtk("mode", modes[2])

plot_mode(modes[1])
#plot_mode(modes[2])
modes

Perturbations#

Let’s add a minor perturbation and calculate the effective refractive index using perturbation theory:

epsilons_p = ["core" => -0.1im, "box" => -0.0im, "slab" => -0.0im, "clad" => -0.0im]
ε_p(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons_p)[tag]

println(perturbed_neff(modes[1], ε_p ∘ τ))

For comparison, we also calculate directly the effective refractive index:

modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1)
println(n_eff(modes_p[1]))
plot_mode(modes_p[1])
Contents