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
Show 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)
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)
Show code cell output
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5900 nodes
Info : 11829 elements
Info : Done reading 'mesh.msh'
2.390356644284898 + 0.0im
2.407726361006542
- 0.00022984678897567025im
-3.638561559087431
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5923 nodes
Info : 11875 elements
Info : Done reading 'mesh.msh'
2.3903569330278445 + 0.0im
2.3977919904522276
- 7.875856704553493e-6im
-5.103702194132172
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5925 nodes
Info : 11879 elements
Info : Done reading 'mesh.msh'
2.3903567025302914 + 0.0im
2.394487593083526
- 3.2575447934408375e-7im
-6.487109603731899
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5911 nodes
Info : 11851 elements
Info : Done reading 'mesh.msh'
2.3903563314332255 + 0.0im
2.3929874345131394
- 1.4781411723613153e-8im
-7.830284085934665
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5914 nodes
Info : 11857 elements
Info : Done reading 'mesh.msh'
2.3903563311984257 + 0.0im
2.3921795616166484
- 7.121592822023824e-10im
-9.147422860794228
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5911 nodes
Info : 11851 elements
Info : Done reading 'mesh.msh'
2.3903563314332614 + 0.0im
2.3916945897787216
- 3.58690131110657e-11im
-10.445280572203725
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5892 nodes
Info : 11813 elements
Info : Done reading 'mesh.msh'
2.3903562047159586 + 0.0im
2.391380481145199
- 1.8918919592036673e-12im
-11.72310366860094
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5908 nodes
Info : 11845 elements
Info : Done reading 'mesh.msh'
2.3903566734068202 + 0.0im
2.3911659352634143
- 9.949501624514985e-14im
-13.002198672736064
Info : Reading 'mesh.msh'...
Info : 25 entities
Info : 5904 nodes
Info : 11837 elements
Info : Done reading 'mesh.msh'
2.3903566757325736 + 0.0im
2.391012362185479
- 4.8889902285710634e-15im
-14.310780830742281
9-element Vector{ComplexF64}:
2.407726361006542 - 0.00022984678897567025im
2.3977919904522276 - 7.875856704553493e-6im
2.394487593083526 - 3.2575447934408375e-7im
2.3929874345131394 - 1.4781411723613153e-8im
2.3921795616166484 - 7.121592822023824e-10im
2.3916945897787216 - 3.58690131110657e-11im
2.391380481145199 - 1.8918919592036673e-12im
2.3911659352634143 - 9.949501624514985e-14im
2.391012362185479 - 4.8889902285710634e-15im
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]
Show code cell output
5-element Vector{Float64}:
-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)
