diff --git a/docs/julia/waveguide bent.jl b/docs/julia/waveguide bent.jl new file mode 100644 index 00000000..2a38bf04 --- /dev/null +++ b/docs/julia/waveguide bent.jl @@ -0,0 +1,113 @@ +# --- +# jupyter: +# jupytext: +# custom_cell_magics: kql +# formats: jl:percent,ipynb +# text_representation: +# extension: .jl +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.11.2 +# kernelspec: +# display_name: base +# language: julia +# name: julia-1.9 +# --- + +# %% [markdown] +# # Mode solving + +# %% [markdown] +# ```{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) +# ``` + +# %% tags=["hide-output", "thebe-init"] +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 + +radius = 100000 +wg_width = 0.5 +wg_thickness = 0.22 +sim_width = 3 +sim_height = 4 +core = shapely.geometry.box(radius - wg_width / 2, 0, radius + wg_width / 2, wg_thickness) + +env = shapely.geometry.box( + radius - sim_width / 2, + -sim_height / 2, + radius + sim_width / 2, + sim_height / 2, +) + +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.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", +) + +# %% tags=["remove-stderr", "hide-output", "thebe-init"] +using Gridap +using Gridap.Geometry +using Gridap.Visualization +using Gridap.ReferenceFEs +using GridapGmsh +using GridapMakie, CairoMakie + +using Femwell.Maxwell.Waveguide + +CairoMakie.inline!(true) + +# %% tags=["remove-stderr"] +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, "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 = 2, order = 1, radius = radius) +println(n_eff(modes[2])) +# write_mode_to_vtk("mode", modes[2]) + +plot_mode(modes[2], absolute = true) +#plot_mode(modes[2]) +modes + +# %% [markdown] +# For comparison, we the effective refractive index of a straight waveguide: + +# %% tags=["remove-stderr"] + +modes_p = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1) +println(n_eff(modes_p[2])) +plot_mode(modes_p[2], absolute = true) + +# %% diff --git a/docs/julia/waveguide.jl b/docs/julia/waveguide.jl index 07234f1b..9e80f3ce 100644 --- a/docs/julia/waveguide.jl +++ b/docs/julia/waveguide.jl @@ -32,18 +32,12 @@ clip_by_rect = pyimport("shapely.ops").clip_by_rect OrderedDict = pyimport("collections").OrderedDict mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict -radius = 5 wg_width = 0.5 wg_thickness = 0.22 slab_width = 3 slab_thickness = 0.11 -core = shapely.geometry.box(radius - wg_width / 2, 0, radius + wg_width / 2, wg_thickness) -slab = shapely.geometry.box( - radius - slab_width / 2, - 0, - radius + slab_width / 2, - slab_thickness, -) +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( @@ -87,16 +81,16 @@ model = GmshDiscreteModel("mesh.msh") labels = get_face_labeling(model) -epsilons = ["core" => 3.5^2 + 0.001im, "slab" => 3.5^2, "box" => 1.444^2, "clad" => 1.0^2] +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, radius = 5) +modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 1, order = 1) println(n_eff(modes[1])) -write_mode_to_vtk("mode", modes[1]) +# write_mode_to_vtk("mode", modes[2]) plot_mode(modes[1]) #plot_mode(modes[2]) @@ -117,7 +111,8 @@ println(perturbed_neff(modes[1], ε_p ∘ τ)) # %% tags=["remove-stderr"] -#modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1) -#println(n_eff(modes_p[1])) +modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1) +println(n_eff(modes_p[1])) +plot_mode(modes_p[1]) # %% diff --git a/src/Maxwell/Waveguide/Waveguide.jl b/src/Maxwell/Waveguide/Waveguide.jl index a36ccf8e..70bffb84 100644 --- a/src/Maxwell/Waveguide/Waveguide.jl +++ b/src/Maxwell/Waveguide/Waveguide.jl @@ -99,46 +99,58 @@ function calculate_modes( Ω = Triangulation(model) dΩ = Measure(Ω, 2 * order) - μ_r = 1 - radius_factor(x) = (1 + x[1] / radius)^2 - twothree = TensorValue([1 0 0; 0 1 0]) - lhs((u1, u2), (v1, v2)) = - ∫( - 1 / μ_r * (curl(u1) ⊙ curl(v1) / k0^2 + ∇(u2) ⊙ v1) + - radius_factor * ( - -u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ v1 + - u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ ∇(v2) - - u2 * VectorValue(0, 0, 1) ⋅ ε ⋅ VectorValue(0, 0, 1) * v2 * k0^2 - ), - )dΩ - rhs((u1, u2), (v1, v2)) = ∫(-1 / μ_r * u1 ⊙ v1 / k0^2)dΩ - - radius_function(x) = x[1] - r = interpolate_everywhere(radius_function, V2) - lhs_radial((u1, u2), (v1, v2)) = - ∫( - -(r / radius * curl(u1) ⋅ curl(v1)) + - (radius / r * gradient(u2) ⋅ v1) + - (k0^2 * ε * r / radius * u1 ⋅ v1) - (radius / r * gradient(u2) ⋅ gradient(v2)) + - (k0^2 * ε * radius / r * u2 * v2), - )dΩ - - rhs_radial((u1, u2), (v1, v2)) = - ∫((radius / r * u1 ⋅ v1) - (radius / r * u1 ⋅ gradient(v2)))dΩ - epsilons = ε(get_cell_points(Measure(Ω, 1))) k0_guess = isnothing(k0_guess) ? k0^2 * maximum(maximum.(maximum.(real.(epsilons)))) * 1.1 : k0_guess + lhs, rhs = if radius == Inf + μ_r = 1 + twothree = TensorValue([1 0 0; 0 1 0]) + lhs_straight((u1, u2), (v1, v2)) = + ∫( + 1 / μ_r * (curl(u1) ⊙ curl(v1) / k0^2 + ∇(u2) ⊙ v1) + ( + -u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ v1 + + u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ ∇(v2) - + u2 * VectorValue(0, 0, 1) ⋅ ε ⋅ VectorValue(0, 0, 1) * v2 * k0^2 + ), + )dΩ + rhs_straight((u1, u2), (v1, v2)) = ∫(-1 / μ_r * u1 ⊙ v1 / k0^2)dΩ + + lhs_straight, rhs_straight + else + radius_function(x) = x[1] + r = interpolate_everywhere(radius_function, V2) + lhs_radial((u1, u2), (v1, v2)) = + ∫( + -(r / radius * curl(u1) ⋅ curl(v1)) + + (radius / r * gradient(u2) ⋅ v1) + + (k0^2 * ε * r / radius * u1 ⋅ v1) - + (radius / r * gradient(u2) ⋅ gradient(v2)) + + (k0^2 * ε * radius / r * u2 * v2), + )dΩ + + rhs_radial((u1, u2), (v1, v2)) = + ∫((radius / r * u1 ⋅ v1) - (radius / r * u1 ⋅ gradient(v2)))dΩ + + lhs_radial, rhs_radial + end + + assem = Gridap.FESpaces.SparseMatrixAssembler(U, V) - A = assemble_matrix(lhs_radial, assem, U, V) - B = assemble_matrix(rhs_radial, assem, U, V) + A = assemble_matrix(lhs, assem, U, V) + B = assemble_matrix(rhs, assem, U, V) if all(imag(A.nzval) .== 0) && all(imag(B.nzval) .== 0) A, B = real(A), real(B) end vals, vecs = eigs(A, B, sigma = k0_guess, nev = num) - vecs[num_free_dofs(V1)+1:end, :] ./= 1im * sqrt.(vals)' / k0^2 + + if radius == Inf + vecs[num_free_dofs(V1)+1:end, :] ./= 1im * sqrt.(vals)' / k0^2 + else + vecs[num_free_dofs(V1)+1:end, :] ./= 1im * sqrt.(vals)' + end + return [ Mode( @@ -157,7 +169,13 @@ function plot_field(field) display(fig) end -function plot_mode(mode::Mode; vertical = false, vectors = false, same_range = false) +function plot_mode( + mode::Mode; + vertical = false, + vectors = false, + same_range = false, + absolute = false, +) Ω = get_triangulation(mode.E) model = get_active_model(Ω) labels = get_face_labeling(model) @@ -190,8 +208,8 @@ function plot_mode(mode::Mode; vertical = false, vectors = false, same_range = f plt = plot!( ax, Ω, - real((E(mode) ⋅ vector)), - colorrange = (-colorrange, colorrange), + (absolute ? abs : real)((E(mode) ⋅ vector)), + #colorrange = (-colorrange, colorrange), ) wireframe!(fig[x, y], ∂Ω, color = :black) Colorbar(fig[x+!vertical, y+vertical], plt, vertical = vertical)