Merge pull request #100 from HelgeGehring/cylindrical_mode_solver

Cylindrical mode solver

docs/_toc.yml

@@ -43,6 +43,7 @@ parts:
   - caption: Julia Examples
     chapters:
     - file: julia/waveguide.jl
+    - file: julia/waveguide_bent.jl
     - file: julia/waveguide_anisotropic.jl
     - file: julia/waveguide_overlap_integral.jl
     - file: julia/lithium_niobate_phase_shifter.jl

docs/julia/waveguide bent.jl

@@ -0,0 +1,131 @@
+# ---
+# 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 for bent waveguides
+
+# %% [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 = 1
+wg_width = 0.5
+wg_thickness = 0.22
+sim_left = 0.5
+sim_right = 3
+sim_top = 1
+sim_bottom = 3
+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.015, "distance" => 0.5),
+    "slab" => Dict("resolution" => 0.015, "distance" => 0.5),
+)
+
+mesh = mesh_from_OrderedDict(
+    polygons,
+    resolutions,
+    default_resolution_max = 0.2,
+    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.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]
+
+
+#dΩ = Measure(Ω, 1)
+τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
+pml_x = x -> 0
+pml_y = x -> 0
+pml_x = x -> 0.1 * max(0, x[1] - (radius + wg_width / 2 + sim_right))^2
+pml_y = x -> 0.1 * max(0, -x[2] - sim_bottom)^2
+
+modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1)
+modes = calculate_modes(
+    model,
+    ε ∘ τ,
+    λ = 1.55,
+    num = 2,
+    order = 1,
+    radius = radius,
+    pml = [pml_x, pml_y],
+    k0_guess = modes[1].k,
+)
+println(n_eff(modes[1]))
+println(log10(abs(imag(n_eff(modes[1])))))
+# write_mode_to_vtk("mode", modes[2])
+
+plot_mode(modes[1], 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[1]))
+#plot_mode(modes_p[1], absolute = true)
+
+# %%

docs/julia/waveguide.jl

@@ -32,23 +32,30 @@ clip_by_rect = pyimport("shapely.ops").clip_by_rect
 OrderedDict = pyimport("collections").OrderedDict
 mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict
 
-wg_width = 2.5
-wg_thickness = 0.3
-core = shapely.geometry.box(-wg_width / 2, 0, +wg_width / 2, wg_thickness)
-env = shapely.affinity.scale(core.buffer(5, resolution = 8), xfact = 0.5)
+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))
+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 = 10,
+    default_resolution_max = 1,
     filename = "mesh.msh",
 )
 
@@ -74,27 +81,27 @@ model = GmshDiscreteModel("mesh.msh")
 
 labels = get_face_labeling(model)
 
-epsilons = ["core" => 1.9963^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 = 2, order = 1)
+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])
+#plot_mode(modes[2])
 modes
 
 # %% [markdown]
 # ## Perturbations
 # Let's add a minor perturbation and calculate the effective refractive index using perturbation theory:
 
 # %% tags=["remove-stderr"]
-epsilons_p = ["core" => -0.1im, "box" => -0.0im, "clad" => -0.0im]
+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 ∘ τ))
@@ -106,5 +113,6 @@ println(perturbed_neff(modes[1], ε_p ∘ τ))
 
 modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1)
 println(n_eff(modes_p[1]))
+plot_mode(modes_p[1])
 
 # %%

src/Maxwell/Waveguide/Waveguide.jl

@@ -74,6 +74,7 @@ function calculate_modes(
     radius::Real = Inf,
     k0_guess::Union{Number,Nothing} = nothing,
     metallic_boundaries = String[],
+    pml::Union{Vector{Function},Nothing} = nothing,
 )
     if count(isnothing, [k0, λ]) != 1
         throw(ArgumentError("Exactly one of k0,λ must be defined"))
@@ -99,24 +100,59 @@ 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Ω
-
     epsilons = ε(get_cell_points(Measure(Ω, 1)))
     k0_guess =
         isnothing(k0_guess) ? k0^2 * maximum(maximum.(maximum.(real.(epsilons)))) * 1.1 :
-        k0_guess
+        k0_guess^2
+
+    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
+        V_pml = TestFESpace(model, ReferenceFE(lagrangian, Float64, 2))
+        pml_f = interpolate(pml[1], V_pml)
+        pml_g = interpolate(pml[2], V_pml)
+
+        s_r = 1 - 1im * gradient(pml_f) ⋅ VectorValue(1, 0)
+        s_y = 1 - 1im * gradient(pml_g) ⋅ VectorValue(0, 1)
+        γ_θ = 1 - 1im * pml_f / (x -> x[1])
+
+        d_θ = s_r * s_y / γ_θ
+        D_t =
+            γ_θ *
+            (s_y / s_r * TensorValue(1, 0, 0, 0) + s_r / s_y * TensorValue(0, 0, 0, 1))
+
+        radius_function(x) = x[1]
+        r = interpolate_everywhere(radius_function, V2)
+        lhs_radial((u1, u2), (v1, v2)) =
+            ∫(
+                -(1 / d_θ * r / radius * curl(u1) ⋅ curl(v1)) +
+                (radius / r / (γ_θ * γ_θ) * (D_t ⋅ gradient(u2)) ⋅ v1) +
+                (k0^2 * ε * r / radius * (D_t ⋅ u1) ⋅ v1) -
+                (radius / r / (γ_θ * γ_θ) * (D_t ⋅ gradient(u2)) ⋅ gradient(v2)) +
+                (k0^2 * ε * d_θ * radius / r * u2 * v2),
+            )dΩ
+
+        rhs_radial((u1, u2), (v1, v2)) =
+            ∫(
+                (radius / r * (D_t ⋅ u1) ⋅ v1) -
+                (radius / r / (γ_θ * γ_θ) * (D_t ⋅ u1) ⋅ gradient(v2)),
+            )dΩ
+
+        lhs_radial, rhs_radial
+    end
+
 
     assem = Gridap.FESpaces.SparseMatrixAssembler(U, V)
     A = assemble_matrix(lhs, assem, U, V)
@@ -125,7 +161,13 @@ function calculate_modes(
         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(
@@ -144,7 +186,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)
@@ -177,8 +225,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)